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Zusammenfassung 


Aufgrund ihrer Kombination aus Festigkeit und Zahigkeit werden poly- 
kristalline Metalle oft fiir Strukturbauteile verwendet, welche anspruchs- 
vollen Betriebsbedingungen, z. B. zyklischen mechanischen Belastungen, 
ausgesetzt sind. Wegen der hohen Sicherheitsvorgaben sowie der ökolo- 
gischen und wirtschaftlichen Anforderungen ist es unerlässlich für eine 
sichere und kosteneffiziente Konstruktion, die eingesetzten Werkstoffe 
für die jeweiligen Belastungsbedingungen zu optimieren. Insbesondere 
bei zyklischer mechanischer Beanspruchung werden die Eigenschaf- 
ten des Werkstoffes stark durch die zugrunde liegende Mikrostruktur 
beeinflusst. Darum muss diese Wechselwirkung bei der Gestaltung 
einer Komponente berücksichtigt werden. Die numerische Homoge- 
nisierung bietet eine Möglichkeit den Einfluss der Mikrostruktur auf 
die makroskopischen mechanischen Eigenschaften von polykristallinen 
Werkstoffen abzubilden. Die vorliegende Arbeit beschäftigt sich mit der 
Erzeugung von Abbildungen polykristalliner Mikrostrukturen welche 
als Eingangsgröße für die numerische Homogenisierung dienen können. 
Des Weiteren wird für die robuste und effiziente Parametrisierung des 
zugrunde liegenden mikromechanischen Modells anhand makroskopi- 
scher Experimente ein Optimierungsansatz vorgeschlagen, welcher auf 
Methoden des maschinellen Lernens beruht. 

Die Abbildung der Mikrostruktur, auch synthetische Mikrostruktur 
genannt, muss relevante Eigenschaften, wie zum Beispiel die Korn- 
größenverteilung und kristallografische Textur des realen Polykristalls 
widerspiegeln. Hierfür werden in dieser Arbeit sogenannte Laguerre- 
Tessellierungen verwendet, welche die Domäne der Mikrostruktur in 
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einzelnen Zellen zerlegt. Die Unterteilung erfolgt auf Basis von Nu- 
klei, welche zuvor in der Domäne zufällig verteilt wurden und den 
Tessellierungsgewichten. Durch das Lösen eines konvexen Optimie- 
rungsproblems können die Gewichte so bestimmt werden, dass die 
Zellen eine vorgeschriebene Größe haben. Hierdurch können syntheti- 
sche Mikrostrukturen erzeugt werden, welche experimentell bestimmte 
Korngrößenverteilungen abbilden. Ebenso können mehrphasige synthe- 
tische Polykristalle mit vorgegebenen Volumenanteilen erzeugt werden. 
Moderne, gradientenbasierte Löser werden im Hinblick auf die Anzahl 
der Iteration und die Gesamtlaufzeit, welche sie zur Lösung des Opti- 
mierungsproblems benötigen, miteinander verglichen. 

Im nächsten Schritt wird eine Methodik entwickelt um den berechneten 
Laguerre-Zellen kristallografische Orientierungen zuzuweisen. Hierzu 
wird eine Kostenfunktion definiert, welche den Unterschied zwischen 
vorgegebenen und aus der synthetischen Mikrostruktur berechneten 
Koeffizienten einer Fourier-Reihenentwicklung der kristallographischen 
Orientierungsverteilungsfunktion beschreibt. Zur Minimierung des Un- 
terschiedes zwischen diesen sogenannten Texturkoeffizienten wird ein 
Gradientenabstiegsverfahren verwendet, welches iterativ die kristallo- 
graphischen Orientierungen der synthetischen Mikrostrukturen anpasst. 
Mit diesem Ansatz wird die Möglichkeit untersucht, isotropes sowie 
anisotropes mechanisches Verhalten zu reproduzieren, sowohl bei zu- 
grundeliegenden gleich- als auch log-normal-verteilten Korngrößen. 
Aufgrund der optimierten Orientierungen reproduziert die synthetische 
Mikrostruktur das gewtinschte mechanische Verhalten mit vergleichs- 
weise wenigen Körnern, d.h. kleinen Volumenelementen. 

Um auf Basis der zuvor generierten synthetischen Mikrostrukturen das 
mechanische Verhalten vorherzusagen, muss das verwendete mikrome- 
chanische Materialmodell parametrisiert werden. Es wird vorgeschlagen 
mithilfe von makroskopischen polykristallinen Experimenten die mikro- 
mechanischen Parameter zu identifizieren. Hieraus ergibt sich ein inver- 
ses Optimierungsproblem, wobei die Auswertung der sich daraus erge- 
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bende Zielfunktion ressourcenintensiv ist. Zur Minimierung der nötigen 
Funktionsaufrufe wird die Bayes’sche Optimierungsstrategie vorgeschla- 
gen, welche ein statistisch informiertes Ersatzmodell der Kostenfunktion 
aufbaut und mithilfe geeigneter Erfassungsfunktionen die Erforschung 
und Ausnutzung ausbalanciert. Die Leistungsfahigkeit und Robustheit 
der Bayes’schen Strategie wird untersucht und mit anderen, für ähnliche 
Probleme verwendete, Optimierungsalgorithmen verglichen. 
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Summary 


Due to their combination of strength and toughness, polycrystalline 
metals are often used for structural components that are subjected to 
demanding operating conditions, e.g. cyclic mechanical loads. Owing 
to the high safety requirements as well as environmental and economic 
restrictions, tailoring these materials to their respective loading condi- 
tions for a safe and cost-efficient design seems imperative. Especially 
in the case of cyclic mechanical loading, the behavior of the material 
is strongly influenced by the underlying microstructure. Thus, this 
interaction must be taken into account when designing a component. 
Computational homogenization offers a way to capture the influence of 
the microstructure on the macroscopic mechanical properties of poly- 
crystalline materials. In this work we address generating representa- 
tions of polycrystalline microstructures which can be used as input for 
computational homogenization. Furthermore, we propose using an 
optimization approach based on machine-learning, for a robust and 
efficient parameterization of the underlying micromechanical model, 
based on macroscopic experiments. 

The representation of the microstructure, also called synthetic microstruc- 
ture, must reflect relevant properties, e.g., grain size distribution and 
crystallographic texture, of the real polycrystal. For this purpose we 
use so-called Laguerre tessellations, which decomposes the synthetic 
microstructure domain into individual cells. This decomposition uses, 
previously in the domain randomly distributed, nuclei and tessellation 
weights. By solving a convex optimization problem, the weights can be 
determined so that the cells have a predefined size. This allows generat- 
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ing synthetic microstructures, that reproduce experimentally observed 
grain size distributions. Similarly, multi-phase synthetic polycrystals 
with given volume fractions can be generated. We compare different 
modern gradient-based solvers in terms of the number of iterations and 
overall run time they need to solve the optimization problem. 

In a next step, we develop a methodology to assign crystallographic 
orientations to the Laguerre cells. For this purpose, we define a cost 
function, describing the difference between given and from the synthetic 
microstructure computed coefficients of a Fourier series expansion of 
the crystallographic orientation distribution function. For minimizing 
the difference between these, so-called texture coefficients, we use a 
gradient descent method, which iteratively adjusts the crystallographic 
orientations of the synthetic microstructures. We investigate the ca- 
pability of the method to reproduce isotropic as well as anistropic 
mechanical behavior, both with underlying uniform and log-normal 
grain size distributions. Due to the optimized orientations, the synthetic 
microstructure reproduces the desired mechanical behavior with few 
grains, i.e., small volume elements. 

To predict the mechanical behavior based on the previously generated 
synthetic microstructures, the micromechanical material model must 
be parameterized. We propose using macroscopic polycrystalline ex- 
periments to identify the micromechanical parameters. This results in 
an inverse optimization problem with an objective function, which is 
resource intensive to evaluate. We propose to use a Bayesian optimiza- 
tion strategy in order to minimize the necessary function evaluations, 
which builds a statistically informed surrogate model of the cost function 
and uses appropriate acquisition functions to balance exploration and 
exploitation. We investigate the performance and robustness of the 
Bayesian optimization approach and compare it to other optimization 
algorithms used for similar problems. 
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Chapter 1 


Introduction 


1.1 Motivation 


Meeting the high safety, environmental and customer demands in the 
automotive industry can be quite challenging. Especially, since the 
loading conditions can be complex and environmental restrictions re- 
quire a light-weight design in order to reduce energy consumption. 
To construct components fulfilling these restrictions, suitable materials 
have to be identified. 

Low-alloyed, high-strength carbon steels are often employed because 
of their high strength and toughness while being cost-effective, see Hor- 
vath (2021) for an overview on applications in the automotive sector. 
Designing components made from these materials in a safe way, while 
taking into account the need for weight reduction, requires considering 
the loading conditions which the component has to face. Approximately 
90% of metallic components fail due to fatigue (Richard et al., 2014), 
which describes failure under cyclic loading conditions below the static 
failure limit. To prevent fatigue upon everyday use, it is necessary to 
assess the cyclic mechanical behavior, including the fatigue limit, of the 
component under consideration. 

Typically, the fatigue limit is determined via mechanical experiments. 
A pertinent engineering tool to assess the fatigue properties of a com- 
ponent are Wöhler diagrams, which show the applied strain or stress 
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amplitude over the endured cycles. For the same design and material of 
a component, there is a stochastic influence on lifetimes, i.e., for the same 
loading conditions, the cycles to failure scatter in the Wöhler diagram. 
The cause for this scattering lies in the inhomogeneous microstructure 
of a component’s material, which dictates the deformation behavior 
under given loading conditions, see Sec. 2.2 for a discussion about 
polycrystalline metals and Francois (1989) as well as Krupp (2007) for 
further insights. 

Thus, designing metallic components to endure cyclic loading conditions 
requires taking the underlying microstructure into account. Fig. 1.1 
shows an image of a ferritic microstructure (Arnaudov et al., 2020), 
obtained by Electron Back Scatter Diffraction (EBSD) methods. 


Figure 1.1: Image of a polycrystalline microstructure obtained by Electron Backscatter 
Diffraction (EBSD) (Arnaudov et al., 2020) and postprocessed by MTex (Nolze and 
Hielscher, 2016). 


Each crystallite is colored according to its orientation with respect to 
the laboratory system. We observe, even for a single phase polycrys- 
tal, a spatial distribution of crystallographic orientations within the 
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microstructure. As the mechanical properties depend on the crystallite 
orientation, there is a spatial distribution of mechanical properties within 
the polycrystal, i.e., a spatial variation of the deformation behavior. 
This heterogeneity in the microstructure leads to the scatter typically 
observed in fatigue tests (Krupp, 2007). 

To reduce time- and cost-intensive experiments, required to capture the 
stochastic influence on fatigue caused by the microstructure, analytical 
and computational homogenization techniques may be used. These 
methods allow determining the effective mechanical behavior based 
on available microstructure data. Whereas being computational more 
demanding than analytical approaches, numerical full-field homoge- 
nization schemes provide access to local fields, e.g., stresses and strains 
as well as plastic deformations, forming in the microstructure. These 
fields serve as input to predict the fatigue behavior of a material under 
given loading (Schäfer et al., 2019a;b; Natkowski et al., 2021a;b). 

The above methods require input data in terms of a computational 
domain, representing the microstructure, as well as a material model, 
describing the deformation behavior on the microscale. This thesis 
provides novel tools to supply this input for polycrystalline materi- 
als to computational homogenization software, i.e., methods to create 
computational cells and identify parameters for the micromechanical 
material model. 
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1.2 Outline 


In Chap. 2 we introduce the foundations of this work. Starting from the 
basics of continuum mechanics, we introduce a pertinent material model 
to capture the described deformation behavior. Additionally we provide 
an overview of the basics of homogenization and the necessary input 
for computational homogenization. Throughout Chap. 3 to Chap. 5 we 
introduce the novelties forming the basis of this work. 

In Chap. 3 we build upon the work by Bourne et al. (2020). They enabled 
the fast computation of Laguerre tessellations with cells of prescribed 
volume fraction by solving a convex optimization problem. We harness 
the potential of modern gradient-type-solvers to speed up solving this 
convex problem and compare them for a variety of applications. In 
addition, we use Anderson acceleration to reduce the computation 
time when generating so-called centroidal tessellations, i.e., tessellations 
where the seed of a cell coincides with its centroid. The combination 
of modern solvers and Anderson acceleration enables the generation 
of (centroidal) polycrystalline microstructures of industrial complexity 
within seconds up to a few minutes, improving upon other reported 
polycrystalline microstructure generators. Indeed, we show that with 
this method it is possible to generate a polycrystalline morphology which 
matches experimentally observed statistics, e.g., grain size distribution, 
with little computational effort. 

Extending the previously proposed method to generate tessellations, rep- 
resenting polycrystalline microstructures, we propose a novel method 
to equip the resulting cells with suitable crystallographic orientations in 
Chap. 4. Using the texture coefficients, i.e., the tensorial coefficients 
of a Fourier series expansion of the crystallite orientation distribution 
function, we formulate a cost function which computes the difference 
between prescribed and realized texture coefficients. To minimize this 
difference, we use a gradient descent scheme, iteratively adjusting the 
crystallographic orientations within the synthetic microstructure. We 
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compare the proposed "Texture coefficient Optimization for Prescribing 
orientations” (TOP) method to state-of-the-art approaches. To do so, we 
compare the effective linear-elastic and the non-linear plastic behavior 
for a number of realizations and a varying number of grains as well as 
for the uniform and a textured case. Last but not least we investigate the 
capability of the proposed method to reproduce a desired mechanical 
behavior with an additional underlying grain size distribution. 

Chap. 5 applies the Bayesian optimization approach, typically used in 
machine learning, to the problem of identifying the crystal plasticity 
parameters based on polycrystalline experiments. Because extending 
the homogenization software by automatic differentiation strategies 
is difficult (or close to impossible when using commercial software), 
solving this inverse optimization problem requires minimizing a black- 
box function. The Bayesian optimization framework iteratively builds 
up a statistically informed surrogate model of the underlying function 
and uses this information to drive the optimization procedure. Test- 
ing this approach on different search spaces, we gain insight into the 
performance of Bayesian optimization and the energy landscape of the 
underlying black-box function. By comparing the Bayesian approach to 
commonly used algorithms in materials science as well as performant 
derivative-free algorithms, we show the efficiency and robustness of the 
Bayesian optimization when identifying crystal plasticity parameters. 
The methods and main results are summarized in Chap. 6 with a brief 
overview of possible applications. As a closing point, future potential 
fields of research are discussed. 


Chapter 2 


Fundamental concepts 


2.1 Basics of continuum mechanics 


2.1.1 Objectives 


To predict the cyclic behavior of metallic components we need a suitable 
tool set, enabling the description of deformation under arbitrary loading 
conditions. The framework of continuum mechanics provides exactly 
these tools, allowing the rigorous description of body movement, defor- 
mation behavior and their thermodynamic basis. This section is based 
on Sect. 1 and 2 of (Šilhavý, 1997; Haupt, 2013), Sect. 2 and 3 in (Bertram, 
2012) and the work by Truesdell and Noll (2004), without any claim of 
originality or completeness. 


2.1.2 Kinematics 


Kinematics refers to the smooth movement of a material body BC R? in 
time through a d-dimensional Euclidean space (in this thesis we consider 
d = 3). We parameterize the current position of a body 6; at time t by 
the function x : Bo x [0, T] — R4 


x= x(X,t), (2.1) 
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where X € R’ denotes an arbitrary material point position in reference 


configuration Bo. Using the spatial inverse mapping x~! 


X =x} (x,t) (2.2) 


allows identifying the reference position X based on the current place- 
ment. The displacement field u : Bo x [0, T] + R?, defined by 


u(X,t)=x(X,t)- X (2.3) 


describes the difference between the current and the reference placement, 
see Fig. 2.1. 


Figure 2.1: Motion and displacement of a continuum body B. 


Describing a field quantity A in terms of the current placement, denoted 
by Ar(z,t), is called the Eulerian description. Using the reference 
configuration to describe the same quantity, i.e., Ar(X,t), is known 
as the Lagrangian description. 
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These are connected through (Haupt, 2013) 


An(z,t) =Arlx "(z,t),t), (2.4) 
Ar(X, t) = Ag(x(X, t),t). (2.5) 


Following Truesdell and Noll (2004), the derivative of A with respect to 
time t yields, for the Lagrangian description, 


Ar(X,t) = za t), (2.6) 


and, in the Eulerian description, 


: _ OAR OAR , 
Ap(z,t) = pp © F By OP 2, (2.7) 
where t = v(x,t) denotes the velocity of the material point. For 


the reader’s convenience we will drop the subscripts E and L in 
the following and indicate the description by the arguments (x,t) 
and (X,t), respectively. 
The deformation gradient F : Bo x [0, T] + Rd, defined by 
F(X,t) = OX (x t) (2.8) 
, Ax rt) 
is a tensorial field quantity, which describes the variation of deforma- 
tion in the neighborhood of a material point from reference to current 
configuration. Thus, it maps the infinitesimal line elements dX in the 
reference placement By to the corresponding line elements dx in the 
current configuration B; 
de=F-dX. (2.9) 


2 Fundamental concepts 


Analogous explicit formulas for the conversion of both area da and 
volume dv elements undergoing the deformation are available 


da = det(F)F~7 . dA (2.10) 
dv = det(F)dV. (2.11) 


For a rigid body movement in three dimensions without rotation, the 
corresponding deformation gradient computes to F = Id, where Id 
denotes the identity on rank two tensors. 

We quantify the deformation in an infinitesimal neighborhood around 
a material point X, e.g., the change of length and angles between line 
elements, by the Green’s strain (Haupt, 2013) 


EF = - (FTF — 1d). (2.12) 


NI = 


Using the displacement gradient H : Bo x [0, T] > R*@, defined by 


Ou 
A(X,t) = —-(X,t 2.1 
KIN), (2.13) 
Green’s strain reads 
E? = (H+H + HTH). (2.14) 


If the displacements are small w.r.t. the typical dimensions in the 
continuum, i.e., 
Bı ~ Bo (2.15) 


the material and spatial coordinates x and X are approximately equal, 
ie, x & X. Thus, it is not necessary to distinguish between the La- 
grangian and Eulerian description, i.e., 


A(X, t) = A(x, t). (2.16) 
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For this case, the Frobenius norm ||-|| for all H = H(X, t) is 
small (Bertram, 2012) 


\|H|| = \/tr(HHT) <1. (2.17) 


Linearizing the Green’s strain around the rest state H = 0 gives rise to 
the infinitesimal strain tensor (Haupt, 2013) 


E= KE H): (2.18) 


2.1.3 Balance equations 


With suitable tools to describe the motion of a continuous body at 
hand, we turn our attention to the physical principles governing the 
(thermo)mechanical behavior of continuum bodies. These principles 
are formulated in the form of balance equations of the (tensor) quanti- 
ties under consideration. To give an overview of the most relevant 
equations, we formulate the general notation of a balance equation and 
subsequently use this framework to provide specific formulations for the 
relevant fields. For the case of singular surfaces, we refer the interested 
reader to the work by Šilhavý (1997). We assume in the following that 
there are no singular surfaces in the material body. 

It is assumed that the change of a field quantity A is the sum of a 
production term pı and a supply term s, of A in the volume V of 
the body B, and the flux q4 of A across the boundary 0B, see Bertram 
(2012); Truesdell and Noll (2004) for more details. This relation can be 
expressed in integral form (Truesdell and Noll, 2004) 


2 Adv= / (pa + sa) du +f qa: nda. (2.19) 
dt Jg, B, öB, 


Here qa, pa and s, denote tensor fields. pa and s, have the identical 
rank as A, whereas, for gq, the rank is increased by one. Using 
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the divergence theorem in combination with Reynold’s transport 
theorem enables to transform the integral form into a local form 
in regular points (Bertram, 2012) 

On 

= + div(Av) = div(qa) + pa + Sa. (2.20) 
Mass balance The balance of mass emerges in local form by consider- 
ing the density p : B; x [0,7] — R. Thus, the balance of mass reads 


b+ pdiv(v) = 0, (2.21) 


as the production pp, supply s, and flux q, are zero for closed systems. 
For small strain deformations, it is typically assumed that the relation 
p = po holds, i.e., the balance of mass is fulfilled for trivial reasons. 


Linear and angular momentum For the linear momentum density pv, 
the supply term is given by the volume force density b : B; x [0, T] > R“. 
The flux consists of the Cauchy stress tensor o : B; x [0, T] — R*@. Thus, 
the local form of the linear momentum balance reads (Silhavy, 1997) 


pa = div(o) + b, (2.22) 


where a = ä denotes the acceleration, which, in the quasi-static setting 
considered in this work, is zero. The stress tensor o describes the stress 
state in any material point. The symmetry of the stress tensor 


o=o! (2.23) 


is equivalent to the balance of angular momentum (Haupt, 2013). 


Energy and entropy The sum of internal energy density e : By x 
[0,7] — R and kinetic energy density 1/2pv - v is referred to as the 
total energy density. Thus, with the supply of internal heat sources 
w : B, x [0,7] > R, the power of volume forces b - v, the negative 
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heat flux —q and the mechanical power v - ø - n, the second law of 
thermodynamics, i.e., the balance of total energy, reads in the local form 


e+ ag -v) =b-v+w+ div(o? - v) — div(q). (2.24) 


The balance of internal energy results from subtracting the balance 
of linear momentum (2.22) multiplied by the velocity v from equa- 
tion (2.24) (Truesdell and Noll, 2004) 


e=-div()+w+o:e. (2.25) 


We use the entropy density s : B, x [0,7] — R to formulate the local 
form of the entropy balance (Haupt, 2013) 


s= div(qs) + Ps + Ss, (2.26) 


where qs, ps and ss denote the flux, production and supply of entropy, 
respectively. According to the second law of thermodynamics, the 
entropy production is non-negative, i.e., 


Ps > 0. (2.27) 


This conditions restricts the directions of thermodynamic processes (Cole- 
man and Noll, 1974). 
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2.2 Deformation behavior of 
polycrystalline metals 


2.2.1 Crystal structure and deformation mechanisms 


With suitable tools to describe the motion of a continuum at hand, 
we focus on the material which is of primary interest in this work, 
i.e., the class of polycrystalline metals. They can be interpreted as an 
agglomerate of crystals or grains, see Fig.1.1. Each one of these crystals 
has a highly ordered and repetitive atomic lattice (Gottstein, 2013). 
Categorizing the different possible periodic arrangements of crystals 
results in the 14 different Bravais lattice types (Bravais, 1850). Due to its 
importance for industrial applications, we focus on the body-centered 
cubic (bcc) lattice, whose underlying unit cell is a cubic system with 
eight corner atoms and an additional center atom, i.e., in total two atoms 
per unit cell. 

If a load is applied to the crystal lattice, at first it deforms elastically, i.e., 
the deformation is reversible, see Fig. 2.2b for a schematic drawing of an 
elastic deformation. When the load vanishes, the crystal lattice returns 
to its initial configuration, see Fig. 2.2a. It is experimentally observed 
that, if the applied load exceeds a material-specific limit, the structure 
deforms plastically, i.e., the deformation is irreversible. This process is 
depicted in Fig. 2.2, where the applied shear stress r exceeds the limit 
Too resulting in a shift of atoms along a so-called glide plane by one 
inter-atomic distance. 
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jii 


(a) (b) 


AY BH 


(c) (d) 


Figure 2.2: Plastic deformation as the result of shifting the whole upper part at once by 
one atom. (a) T = 0, (b) T < Too, (C) T > Too, (d) T = 0. 
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Using the analogy of a crystal consisting of bubbles, Cottrell (1954) 
estimated the theoretical barrier for this deformation to be 


To = Ë (2.28) 


where u denotes the material’s shear modulus. According to Hull and 
Bacon (2001) this theoretical limit for plastic deformation is, in general, 
several orders of magnitude higher than the experimentally observed 
yield strength. 

An explanation for this striking difference was simultaneously proposed 
by Orowan (1934), Polanyi (1934) and Taylor (1934). Based on the pres- 
ence of one-dimensional defects, so-called dislocations, their movement 
in the lattice, see Fig. 2.3a, is interpreted as the cause of plastic deforma- 
tion. By breaking the atomic bonds of the inserted half plane with one 
of the neighboring planes and bonding again with the following edge, 
see Fig. 2.3b, the dislocation can move in the lattice. This movement 
facilitates the continuous plastic deformation of the crystal lattice and 
explains the reduced resistance against plastic deformation compared to 
equation (2.28). Peierls (1940) and Nabarro (1947) give a discussion of 
the resistance of a crystal lattice against dislocation movement. 
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(a) (b) (c) 


Figure 2.3: Plastic deformation as a result of dislocation movement. (a) Initial configuration, 
(b) Moving dislocation, (c) Final configuration. 


Dislocations generally do not move on planes with arbitrary normal, but 
on predefined systems depending on the crystal lattice (Hull and Bacon, 
2001). The accumulation of dislocation movement on these systems is 
called slip. Those systems which permit dislocation movement are called 
slip systems (Hull and Bacon, 2001). For the body-centered cubic crystal 
there exist 12 primarily activated slip systems. Dislocation movement 
on other systems is only observed in some materials under certain 
temperature conditions (Courtney, 2005), which are not considered in 
this work. Within the crystal, a slip system n is characterized by two 
orthogonal vectors, the slip plane direction m” and the normal n”. 

For uniaxial tension o„ applied to a crystal with slip systems 7, the 
driving force for dislocation movement, i.e., the shear stress, on the slip 
system follows Schmid’s law (Gottstein, 2013) 


T” = On coso” cosà”, (2.29) 


where ¢" and A” denote the angles between the direction of applied 
stress and the slip plane normal and the slip direction, respectively. 
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For uniaxial tension, the shear stress becomes maximal for a slip system 
n with 67 = A" = 45°. 


2.2.2 Single crystal plasticity’ 


To use the powerful computational homogenization methods outlined in 
Sec. 2.3, it is necessary to encapsulate the knowledge of the deformation 
behavior, see Sec. 2.2.1, into a suitable material model. In the case of 
polycrystalline materials, the crystal plasticity model offers a dedicated 
framework to describe the deformation behavior of polycrystalline 
materials accounting for the movements of dislocations. Based on the 
pioneering work by Asaro and Rice (1977), Asaro and Needleman (1985) 
and Peirce et al. (1983), many crystal plasticity models were developed, 
see Roters et al. (2010) for an overview. 

Phenomenological crystal plasticity models use a small number of pa- 
rameters to describe the observed deformation behavior of the crystal- 
lites. Nevertheless, they proved to be able to provide results which, for 
relevant features, agree well with experimental observations (Schafer 
et al., 2019b; Natkowski et al., 2021b). For fatigue behavior, the defor- 
mations are typically small. Thus, we use a crystal plasticity model at 
small-strains, closely following Schäfer et al. (2019a). We assume the 
total strain tensor £ to be decomposed additively 


C= eee (2.30) 


into an elastic and a plastic contribution, denoted by £. and e,, respec- 
tively. We refer to (Simo and Hughes, 2006, Cha. 1) for details. The 
Cauchy stress tensor ø, which is symmetric by conservation of angular 
momentum (Haupt, 2013), is assumed to be linear in the elastic strain, 


1 This subsection is based on excerpts from the publications by Kuhn et al. (2021; 2022) 
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i.e., Hooke’s law 
oe=C:.=C:(e-e,), (2.31) 


involving the fourth-order stiffness tensor C, is assumed to hold. For 
elasto-viscoplasticity, the evolution of plastic strain e, is typically gov- 
erned by a flow rule of the form 


Ep = f(o,2), (2.32) 


where z denotes a vector of additional internal variables, see Chap. 2 
in Simo and Hughes (2006). 

In the framework of crystal plasticity, it is assumed that plastic defor- 
mation, i.e., èp # 0, is caused by the movement of dislocations on 
the corresponding crystallographic slip systems. We consider volume- 
preserving slip mechanisms, i.e., conservative glide. An arbitrary slip 
system n is activated when the resolved shear stress 


rl =o-M" (2.33) 


exceeds a critical value r.. Here, M” denotes the symmetrized Schmid 
tensor of the slip system 77 


M” = ; (m! an" +n" 8m”). (2.34) 
Following Bishop (1953) and Hutchinson (1976), the flow rule (2.32) is 


formulated as the superposition of crystallographic slip rates on the 
individual slip systems Ng 


Ns 
èp = X AM", (2.35) 
n=1 


where ¥” denotes the plastic slip rate for the n-th system. To determine 
the amount of plastic slip on a slip system 7, the theory of elasto- 
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viscoplasticity assumes the plastic slip rate 7" to be a function of the 
resolved shear stress. The flow rule is designed in such a way that it 
captures, both, the Bauschinger effect and the ratcheting behavior, which 
are necessary for describing the cyclic behavior of metals accurately. The 
Bauschinger effect concerns a decreased resistance to plastic deformation 
also when a change of loading directions (e.g., from tension to compres- 
sion) occurs. Ratcheting refers to a consistent accumulation of plastic 
slip under stress-driven cyclic loading conditions. Please note that the 
loading levels which lead to ratcheting under repeated loading will not 
induce a plastic flow of the material under static loading conditions, 
in general. For ratcheting and strain-driven loading, the mean stress 
will decrease for increasing number of cycles, see Farooq et al. (2020) 
and Cruzado et al. (2020) for recent discussions of ratcheting and the 
Bauschinger effect, respectively. 

In this work, the flow rule proposed by Hutchinson (1976) 


"= yo sgn(r” — Xp) (2.36) 


is used. The symbol 7,’ denotes the critical resolved shear stress and c 
refers to the stress exponent. To capture both the Bauschinger effect and 
ratcheting behavior, the original flow rule is augmented by the backstress 
X;|, as proposed by Cailletaud (1992). Alternative approaches for model- 
ing these effects are discussed by Harder (2001) with special focus on the 
coupling of the backstress evolution on different slip systems and the 
existence of the backstress tensor. Please note that the vector z of internal 
variables in equation (2.32) collects the different backstresses X). 

To close the model, the flow rule needs to be complemented by an 
evolution equation for the backstress X". The Armstrong-Frederick (AF) 
model (Frederick and Armstrong, 2007) is a nonlinear extension of the 
model by Prager (1949) and involves a recall term. The model describes 
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the evolution of the backstress by the equation 
XT = AV" — BY"), (2.37) 


where A and B are material parameters with dimensions MPa and 
1, respectively. These parameters represent the direct hardening as 
formulated in the Prager model and the recovery modulus present in 
the extended model, respectively. Bari and Hassan (2000) showed that 
the AF model may overestimate ratcheting encountered in experiments. 
Furthermore, they attribute this shortcoming to the constant ratcheting 
associated to the model. 

Extending a formulation by Chaboche, which is based on a decompo- 
sition of the AF model, Ohno and Wang (1993) proposed a kinematic 
hardening model that reads, in the context of crystal plasticity, 


xX = Aà -B Ge i xn a). (2.38) 
b A / B b 
In the Ohno-Wang (OW) model, the dynamic recovery term is mod- 
ified by a power law with exponent M, which controls the degree 
of non-linearity, to improve the prediction of ratcheting and mean 
stress relaxation. 
According to Schäfer et al. (2019a;b) and Natkowski et al. (2021a;b) 
the above model captures relevant aspects of fatigue. Thus, in the 
following it serves as our model for the deformation behavior of 
polycrystalline microstructure. 
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2.3 Homogenization methods 


2.3.1 Mechanical properties of microstructured materials 


The deformation mechanisms on the microscopic scale, e.g., outlined in 
Sec. 2.2.1 for polycrystalline metals, influence the overall macroscopic 
behavior of the material. Thus, within many applied materials, we 
face different scales when investigating mechanical properties. Let 
us consider the component shown in Fig. 2.4 as an example. On the 
component, or macroscopic, scale, the metallic material is assumed to 
be homogeneous. On the microscopic scale, see Fig. 2.4 for a example 
of the components microstructure, we observe that this not the case. 
Owing to the spatial distribution of crystallographic orientations and the 
direction-dependent mechanical properties of crystallites, the microstruc- 
ture is heterogeneous. These microstructural heterogeneities influence the 
macroscopic mechanical properties and can be crucial for failure of the 
component, see Francois (1989) and Krupp (2007) for an overview in the 
context of fatigue failure. 

Deriving macroscopic properties based on information on the microscale 
is called homogenization whereas localization refers to deducing mi- 
crostructural behavior from macroscopic quantities, e.g., strains and 
stresses. To determine the effective properties of a microstructure, 
we have to consider a domain which is characteristic for the whole 
microstructure. This representative volume element (RVE) (Hill, 1952; 
1963) needs to be sufficiently large to capture all relevant statistics 
of the microstructure, while being sufficiently small compared to the 
overall size of the specimen. We show this concept of scale separation 
schematically in Fig. 2.4. The component scale fmacro is much larger 
than the microstructural scale (micro, which is larger than scale of one 
constituent, £constituent, €.g-, the size of one grain in polycrystalline metals. 
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Lconstituent 


Figure 2.4: Schematic illustration of different lengths on the micro and macroscopic scale. 


Concluding that the mechanical behavior of the RVE is representative 
for the whole material allows determining effective macroscopic fields, 
e.g., strains € and stresses 7, from averaging the fluctuating micro stress 
and strain fields, o(x) and ¢() over the RVE volume V, i.e., 


Qı 
II 


jo dV (2.39) 


E 


> I steav: (2.40) 


To keep this section focused on the main aspects in this work, we 
discuss analytical and computational homogenization methods with 
their respective requirements as well as a advantages and disadvantages 
in the following. For more details regarding micromechanics and homog- 
enization methods we refer to Nemat-Nasser and Hori (1993), Torquato 
(2002) and Gross and Seelig (2017). 
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2.3.2 Analytical approaches 


The works of Voigt (1889) and Reuss (1929) describe the first analytical 
approaches to homogenization. Both authors approximate the effective 
elastic moduli of a polycrystal by volume averaging the elastic constants 
of single phases. 

Eshelby (1957) provided an analytical solution for the strain and stress 
fields of an eigenstrained ellipsoidal inclusion in an infinite homoge- 
neous matrix, for which the strains and stresses inside the inclusion 
are homogeneous. Based on this solution, the mean-field theory pro- 
vides many different homogenization approaches. These methods have 
successfully been applied to determine effective properties of fiber 
reinforced composites (Schemmann et al., 2018; Kehrer et al., 2020) and 
polycrystalline materials (Rieger and Böhlke, 2015). For the linear elastic 
effective behavior of polycrystals, the self-consistent method, proposed 
by Hershey (1954) and Kroner (1958), is used frequently. There exist 
many different extensions to apply this method for non-linear behavior 
as well, e.g., (Hill, 1965; Molinari et al., 1987). In their article, Lebensohn 
et al. (2007) give an overview over various material classes they study 
using the self-consistent approach. Applying a variational framework, 
Hashin and Shtrikman developed bounds on the effective properties of 
heterogeneous microstructures (Hashin and Shtrikman, 1962a;b). These 
bounds are used in a variety of applications, e.g., Fernández and Böhlke 
(2019) applied the Hashin-Shtrikman bounds to the effective elastic 
properties of textured polycrystals. We refer the interested reader to 
Chap. 8 in Gross and Seelig (2017) for a more detailed overview over 
analytical homogenization methods. 

Although some restrictions apply to these techniques, e.g., in terms 
of the possible inclusion geometries considered in Eshelbys solution, 
they are capable of predicting the effective mechanical properties of 
various material classes. Additionally, they are efficient both in terms of 
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computing time and memory requirements, making them attractive for 
industrial applications. 


2.3.3 Computational homogenization 


To circumvent the simplifying assumptions of analytical approaches, 
e.g., regarding the inclusion geometry, computational homogenization 
approaches use spatially resolved microstructures as computational 
cells. In addition to the effective behavior, numerically solving the 
so-called cell problem on the resolved microstructure, delivers the full 
fields of stresses and strains, which emerge. Furthermore, it is possible 
to compute the resulting plastic deformations developing within the 
microstructure on an intracrystalline resolution which can be used as 
input, e.g., to predict crack initiation and propagation in polycrystalline 
metals (Schäfer et al., 2019a;b; Natkowski et al., 2021a;b). 

To derive the emerging fields and thus the effective properties of a 
microstructure, the finite element method (FEM) might be used to solve 
the discretized weak form of the governing equations (Matouš et al., 
2017). This method predicts the effective properties, e.g., of compos- 
ites (Nakamura and Suresh, 1993; Kouznetsova et al., 2001; Borbely et al., 
2001), microstructures with pores (Zhuang et al., 2015) and polycrystals 
(Miehe et al., 1999; Shenoy et al., 2008; Vajragupta et al., 2014). 

In this work we use a method based on the work by Moulinec and 
Suquet (1994; 1998) to numerically determine the effective properties of 
polycrystalline microstructures. This approach is based on the reformu- 
lation of the cell problem into a Lippmann-Schwinger equation and uses 
the fast Fourier transforms (FFT) for its computational resolution. Due 
to its low memory footprint, the efficient FFT implementations and the 
straightforward realization, the method became increasingly popular 
over the past decade. This fueled the development of new algorithms, 
see Schneider (2021) for a recent review on supported discretization 
methods and available solution schemes. Applications of this method 
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include, but are not limited to, damage and fracture mechanics (Chen 
et al., 2019; Ernesti et al., 2020), homogenization of piezoelectricity (Bren- 
ner, 2009; Göküzüm et al., 2019), computing the deformation behavior 
of polycrystalline materials (Lebensohn and Rollett, 2020; Rovinelli et al., 
2020) and composites (Müller et al., 2015; Görthofer et al., 2020). 

We rely on FeelMath (Fraunhofer ITWM, 2021), a commercial FFT based 
micromechanics solver, to determine a solution of our cell problem. 
Fig. 2.5 shows one example for the resulting von Mises stress field in a 
polycrystalline microstructure computed by FeelMath. Nevertheless we 
still need to supply computational cells, i.e., RVEs of the considered poly- 
crystalline microstructure. Additionally we need to identify parameters 
for the material model described in Sec. 2.2.2 to describe the deformation 
behavior on a microscopic scale accurately. 


(a) (b) 


Figure 2.5: Example for computational homogenization. (a) Spatially resolved polycrys- 
talline microstructure, (b) Resulting von Mises stress field at 0.1% macroscopic strain. 
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2.4 Prerequisites for computational 
homogenization of polycrystalline metals 


2.4.1 Polycrystalline computational cells for 
micromechanical simulations 


The morphology of polycrystalline microstructures can be determined 
by dedicated experimental methods, like serial sectioning (Spowart 
et al., 2003; Kubis et al., 2004) and FIB-SEM (Bansal et al., 2006; Groeber 
et al., 2006; Zaefferer et al., 2008)). By combining these with Electron 
Backscatter Diffraction (EBSD) methods (Adams and Olson, 1998; Korte 
et al., 2011) the local orientation of each grain can be resolved. Fig. 1.1 
and Fig. 2.4 show two examples of images obtained by EBSD. 

Directly using these images as an input for computational homoge- 
nization might seem like the most natural approach. However, the 
experimental procedures are quite time and cost intensive (Bargmann 
et al., 2018) and, in addtion, determining a representative image for the 
whole microstructure can be challenging. Due to size limitations of one 
image, it might not be large enough to capture statistical fluctuations 
of the relevant microstructure features. It might even distort these. 
Using the (generally) non-periodic experimental images can lead to 
large required computational cells, e.g., in terms of grain number within 
the image, to accurately predict the effective behavior (Kanit et al., 2003; 
Yang et al., 2019; Schneider, 2021). 

As an alternative, synthetically generated computational cells gained 
popularity, see Fig. 2.6 for a polycrystalline example. These synthetic 
microstructures aim at reproducing statistical features of the real mi- 
crostructure, e.g., grain size and aspect ratio distribution (Gillner and 
Miinstermann, 2017; Prasad et al., 2019; Henrich et al., 2020). They 
offer the benefit of direct control over the microstructural parameters, 
yielding the possibility to tailor the microstructure to certain loading 
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conditions and to capture the statistical influence of fluctuations in the 
microstructural descriptors (Tallman et al., 2020; Tran and Wildey, 2021). 


Figure 2.6: Example of a polycrystalline computational cell. 


There exists a variety of methods to generate synthetic microstructures, 
each tailored to the material class under investigation. In the following 
we focus on polycrystalline metals and refer the interested reader to the 
review article by Bargmann et al. (2018) for a general overview. 

One possible approach is simulating the solidification and grain growth 
occurring during the manufacturing process of a metal. Examples 
include techniques based on Vertex methods (Kawasaki et al., 1989; 
Weygand et al., 1999), cellular automata (Janssens, 2010; Saluja et al., 2012; 
Kühlbach et al., 2016) and Monte Carlo methods (Anderson et al., 1989; 
Radakrishnan and Zacharia, 1995). Using the phase field method (Krill 
IH and Chen, 2002; Nestler, 2005; Hoetzer et al., 2016) allows considering 
physical principles when simulating the forming process of metals. This 
comes at the cost of increased computing time and having to identify 
physically reasonable model parameters. Examples for polycrystalline 
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materials using the phase field approach can be found in the work 
by Graf et al. (2021a;b), simulating the formation of martensite from a 
prior austenite phase for an industrial application. 

Describing the grains as geometric objects dispenses with the computa- 
tional demanding simulation of solidification and grain growth. These 
objects form a tessellation of the computational domain, i.e., there are 
no voids and the objects do not overlap. Assuming the grains to be 
convex polyhedrons, Voronoi tessellations (Aurenhammer, 1991) can be 
used to generate representations of polycrystalline microstructures. This 
approach uses a distinct set of points, called seeds, in the computational 
domain to define a cell by all the points lying closer to one seed than to 
any other seed, see Fig. 2.7 for a two dimensional example. 


(a) (b) 


Figure 2.7: Schematic drawing of a Voronoi tessellation process, following 
Aurenhammer (1991). (a) Seeds in the computational domain, (b) Resulting tessellation. 


For seed positions following a Poisson process, a so-called Poisson- 
Voronoi tessellation, the generated microstructures often do not capture 
relevant properties of the real microstructure, e.g., grain size distribution 
or number of neighboring grains (Dobrich et al., 2004; Luther and Könke, 
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2009). Choosing the seeds based on a precomputed sphere packing 
increases the synthetic microstructures realism (Fritzen et al., 2009; Quey 
et al., 2011). Laguerre tessellations extend the Voronoi tessellations by 
one additional degree of freedom per cell. This so-called tessellation 
weight controls the size of each cell, enabling the generation of polycrys- 
talline microstructures with prescribed volume fractions (Bourne et al., 
2020). For generating synthetic microstructures with non-convex geomet- 
ric models, several approaches are available, e.g., fitting experimental 
data (Alpers et al., 2015; Teferra and Rowenhorst, 2018) or optimization 
procedures to match the shape distributions (Groeber and Jackson, 2014; 
Prasad et al., 2019; Henrich et al., 2020). 

The heterogeneity of a polycrystalline microstructure is not only influ- 
enced by the shape of its grains, but especially by the distribution of 
crystallographic orientations (Adams and Olson, 1998; Kocks et al., 2000). 
Thus, the geometric modeling approaches require an additional step, 
the assignment of crystallographic orientations to each grain model. For 
polycrystalline materials it is important to distinguish between two cases: 
The uniformly distributed orientations and the textured case where 
certain orientations are more likely to be encountered than others. The 
former case results in mechanically isotropic effective properties (Krawi- 
etz, 1999; Bertram et al., 2000; Böhlke and Bertram, 2001), whereas the 
elastic and elasto-plastic behavior in the latter case might, in general, be 
anisotropic. Besides randomly assigning orientations, e.g., see (Gillner 
and Münstermann, 2017; Tu et al., 2019), some dedicated algorithms exist. 
Examples include approaches discretizing the one-point correlation 
function, i.e., the crystallite orientation distribution function (CODF), 
which describes the orientation state of a polycrystal (Tóth and Van 
Houtte, 1970; Melchior and Delannay, 2006). Following their respective 
discrete CODF values, crystallographic orientations are selected and 
assigned to cells in the synthetic microstructure (Chunlei et al., 2000; 
Deka et al., 2006). Building upon this discretization, Eisenlohr and Roters 
(2008) use a combination of stochastic and deterministic approaches 
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to increase sampling efficiency. Biswas et al. (2020) and Prasad et al. 
(2019) iteratively adjust the discretization of the CODF to improve the 
sampling results, whereas Liu et al. (2020a) adjust the parameters 
of the microstructure generator. Quey and Renversade (2018) match 
the special case of a uniform orientation distribution using uniformly 
distributed unit quaternions. Approximating the CODF by superposing 
so-called texture components (Helming, 1996; Wassermann and Grewen, 
2013) allows sampling these components and assigning orientations 
accordingly (Hielscher and Schaeben, 2008). 


2.4.2 Determination of single crystal 
plasticity parameters 


To describe the deformation behavior of polycrystalline metals on a mi- 
crostructural level, we use the crystal plasticity model outlined in Sec. 2.2. 
However, this model still needs to be equipped with suitable parameters. 
If the model parameters cannot be derived using physical principles, we 
will need to determine the parameters so that the model output matches 
experimental observations. Thus, we need to minimize a function which 
describes the difference between experiment and simulation. 

For crystalline metals, ultrasonic testing on a single crystal (Hunting- 
ton, 1947), microtensile testing (Gianola and Eberl, 2009) as well as 
experiments using milled micropillars (Cruzado et al., 2015) and nano- 
indentation tests (Chakraborty and Eisenlohr, 2017; Engels et al., 2019; 
Liu et al., 2020b) give insight into the deformation behavior on a micro- 
scopic scale. Comparing the experimental observations with simulations 
on single crystals, in combination with suitable optimization methods, 
allows determining the crystal plasticity parameters (Chakraborty and 
Eisenlohr, 2017; Engels et al., 2019). 

Using data from polycrystalline experiments allows reducing the experi- 
mental effort associated with single crystal tests. Comparing the macro- 
scopic data with the effective behavior, determined from simulations 
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on polycrystalline RVEs, allows calibrating crystal plasticity parameters, 
see Herrera-Solaz et al. (2014). Schafer et al. (2019a) define a cost 
function based on the difference between experimental and simulated 
stress-strain hysteresis, leading to an optimization problem where the 
objective function is expensive to compute and gradient information is 
not available. Because automatic differentiation (Griewank and Walther, 
2008) might be difficult to integrate into existing homogenization tools, 
derivative free approaches are frequently used in this setting, see Rios 
and Sahinidis (2013) for an overview. However, care has to be taken 
when considering industrial applications, as determining suitable pa- 
rameters are a regular task (Schafer et al., 2019a;b; Natkowski et al., 
2021b; Arnaudov et al., 2020). Indeed, as evaluating the cost function 
requires a full-field simulation it is desirable to minimize the number of 
function evaluations. 
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Chapter 3 


Fast methods for computing 
centroidal Laguerre tessellations 
for prescribed volume fractions 
with applications to 
microstructure generation of 
polycrystalline materials ' 


3.1 Introduction 


In a recent article, Bourne et al. (2020) used a convex variational 
principle for generating Laguerre tessellations modeling polycrys- 
talline microstructures with prescribed volume fractions per grain, see 
Fig. 3.1 for an overview. Their insights led to the development of 
fast techniques for generating polycrystalline microstructures which 
are orders of magnitude faster than other approaches in the litera- 
ture, for instance based on derivative-free optimization (Quey and 
Renversade, 2018), statistical methods (Lautensack, 2008) or cross- 
entropy approaches (Petrich et al., 2019), and more accurate than 


! This chapter is based on the publication by Kuhn et al. (2020) whereas minor 
typographical and formatting changes have been made for cohesion of the manuscript 
and the readers convenience 
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heuristic approaches (Lyckegaard et al., 2011). 

In this chapter we improve upon the proposed method by Bourne et al. 
(2020), as they essentially rely upon black-box optimization methods. 
They use the limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) 
method (Nocedal, 1980) with backtracking for solving the convex 
optimization problem. The backtracking is used to avoid generating 
empty cells during the iterative process, and is also used for the 
corresponding Newton-type methods (Kitagawa et al., 2019). 

In section 3.2, we set up our notation for computing periodic Laguerre 
tessellations and the resulting convex optimization problem. We discuss 
pertinent solvers in section 3.3. Of particular interest here are modern 
non-monotone gradient-type solvers, i.e., the methods of Malitsky and 
Mishchenko (2019) and Barzilai and Borwein (1988), whose application 
to optimal-transport problems appears to be novel. Compared to more 
traditional solvers of (Quasi-)Newton type, they combine low cost per 
iteration with non-monotone convergence behavior, which turns out to 
be beneficial for the overall run-time. Indeed, non-monotone methods 
have gained popularity in recent years (Nesterov, 2004), as the user 
is mostly interested in minimizing the time to solution, and ensuring 
a monotone convergence of the algorithm in question is typically of 
secondary importance. 

For computing centroidal Laguerre tessellations, we review the classical 
Lloyd’s algorithm in subsection 3.2.2. Subsequently, we apply Anderson 
acceleration to Lloyd’s algorithm. Anderson acceleration (Anderson, 
1965) is a dedicated technique for improving the convergence behavior 
of general fixed-point iterations. Its mathematical analysis is more recent, 
establishing it as a particular Quasi-Newton method of multi-secant 
type (Fang and Saad, 2009) and establishing the improved convergence 
behavior also rigorously (Toth and Kelley, 2015; Evans et al., 2020). 

In section 3.4, we study the numerical effectiveness of the proposed 
algorithms for selected problems of increasing complexity. We study the 
dependence of the algorithms on the number of grains and the prescribed 
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grain-size distribution, with and without centering. Furthermore, we 


exhibit their capabilities for polycrystalline materials of industrial 
relevance. 


| Input: unit-cell dimensions, number of grains, grain-size distribution (GSD)) 


Vv 
(Determine target volumes of the individual cells according to GSD | 


v 
| Initialize seed positions for Laguerre tessellation | 


v 
| Compute Laguerre weights realizing the desired grain volumes 


by solving a convex optimization problem 
Update seed position 


Centroids 
= seeds? 


no 


yes 


Output: Polyhedral cells of the polycrystalline microstructure | 


Figure 3.1: General flowchart for computing Laguerre-type models of polycrystalline 
microstructures, Bourne et al. (2020). 


3.2 Periodic Laguerre tessellations 
and centering 


3.2.1 Periodic Laguerre tessellations 


Let Y be a rectangular domain in Rt, which we fix as Y = [0, L4) x 
[0, L2) x ... x [0, La) for positive lengths L;. Y shall serve as our unit 
cell of interest, and we furnish it with periodic boundary conditions. We 
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denote the corresponding periodic distance function by d : Y x Y + Rso. 
For a sequence (z1,22,...,2x) € Y™ of distinct points lying inside Y 
and a given set of real weights w € RX, we denote by Dj, Do,... Dx 
the periodic Laguerre tessellation (Aurenhammer, 1991) of Y with seeds 
(%1,@2,...,0K) and weights (wi, w2,..., wg). Each cell D; is a subset of 
Y defined via 


D; = {z eY | d(x, 2;)? — wi < d(x, xj)? =w; for j€{1,2,- ., KP} : 
(3.1) 
The interiors of the Laguerre cells are mutually disjoint, and the Laguerre 
cells cover the original volume Y. Clearly, adding a constant to each 
single weight does not change the associated Laguerre cells. Notice 
that the seed x; does not need to lie in the corresponding Laguerre cell 
D;. Also, the Laguerre cell may be empty. This becomes clear when 
imagining the weight w; as a "price". Loosely speaking, if the price 
is excessive, the turnover, which is related to the volume of the cell, 
will be zero. 
For the special case of identical weights, wı = wa =... = wg, a Laguerre 
tessellation is called a Voronoi tessellation, see Aurenhammer (1991). For 
this special case, the seeds are always contained in the corresponding 
cell. However, the kinds of tessellations which may be generated as 
Voronoi tessellations are rather limited, for instance in terms of the 
achievable volume-fraction distributions, see section 2.1.1 in Quey and 
Renversade (2018). In contrast, Laguerre tessellations may generate any 
normal tessellation with convex cells in R@ for d > 3, see Lautensack 
and Zuyev (2008). In this context, normal means that if the seeds are in 
general position, each k-face of the tessellation arises as the intersection 
of exactly d — k + 1 cells fork = 0,1,...,d—1. 
For this article, Laguerre tessellations serve as simple geometrical models 


? 


for polycrystalline microstructures, e.g., of metallic materials, alloys, 
ceramics or polycrystalline ice. We refer to section 2.1 in Bargmann et al. 
(2018) for a recent overview. 


36 


3.2 Periodic Laguerre tessellations and centering 


For a given sequence (271, %2,...,uK) € Y* of distinct points of lying in- 
side Y and a given set of volume fractions (¢1, ¢2,...,¢x) € RX, each of 
which is non-negative, and which sum to unity, there is a corresponding 
set of weights (w1, w2,..., wx) € R*,s.t. the corresponding Laguerre 
cells realize the prescribed volume fractions, i.e., 


|Di| 


see Aurenhammer et al. (1998), where | D;| denotes the volume of cell D;. 
The proof of this result is interesting in its own right, as it immediately 
gives rise to a computational technique for computing the associated 
weight vector w. As noted earlier, the component-wise addition of a 
constant to the weight vector does not change the associated Laguerre 
cells. Thus, we define 


We = [went 


K 
Bin o} , (3.3) 


i=l 


the (K — 1)-dimensional subspace of RX with vanishing mean. For 
a fixed sequence (21,22,...,tK) € Y* of distinct points lying in- 
side Y and a given set of volume fractions (¢1,¢2,...,¢«) € R#, 
define the function 


g: Wk >R, 


K K 
w => X (Qi: LiLo- La — | Dil) wi + 5 d(x, x4)? dz, (3.4) 
i=1 i=1 7 Di 


where D; denotes the i-th associated Laguerre cell (3.1). Aurenhammer- 
Hoffmann-Aronov (Aurenhammer et al., 1998) have shown that the 
function g (3.4) is a concave function which admits maximizers. Fur- 
thermore, each critical point (which is a maximum by concavity of g) 
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satisfies the first order optimality condition 
0 = Dg(w)[o] 


in terms of the directional derivative of g at w € Wx in direction v € Wx. 
More explicitly, 


K 
Dg(w)[v] = Lyw + sv)| o = > (¢i: LiLo La |Dil) vi, (35) 


i=1 


i.e., any critical point of g satisfies equation (3.2). Furthermore, it can be 
shown that g is even strictly concave on the set 


{w € Wg ||D:i|>0 forall i=1,2,...,K}, (3.6) 


see Kitagawa et al. (2019). It is clear that g is not strictly concave outside 
of the set (3.6). Indeed, suppose that we have w € Wx, s.t. the i-the 
Laguerre cell is empty. Then, the Laguerre cells remain unchanged if we 
increase the i-th weight w;. In particular, the function g is constant along 
the ray 


1 
w+s ei = gbb) , s>(, 


where e; is the i-th unit vector in R*. For completeness, notice that the 
i-th component of the gradient of g on Wx computes as 


see (3.5). Notice that Vg(w) € Wx automatically, i.e., its mean vanishes. 
The concavity of g (or, equivalently, the convexity of —g) permits using 
established convex optimization techniques as in Nocedal and Wright 
(1999) to be applied, see section 3.3 for details. To finish this paragraph, 
we report on a result in Bourne and Roper (2015). The function g (3.4) is 
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even twice differentiable on W;, with second derivative 
D?g(w)[v,v] =v H(w)v, w,v € Wr, 


in terms of the Hessian matrix H(w) € R“** with entries 


Hy; (w) = 2 d(x;, £4) 


(3.8) 
for i # j, where A;; denotes the area of the face between the i-th and 
j-th cell (which is formally set to zero if the cells do not share a face), 
and H;; = - ya Hij fori=1,2,...,K. 

As gis concave and twice differentiable, the Hessian matrix H(w) is 
negative semidefinite. However, H(w) does not have full rank, as all 
rows and columns have zero mean, leading to a zero eigenvalue. In 
particular, it is not invertible. As a countermeasure, we introduce the 


matrix 


~ tr(H(w)) | or 
with the vector 1 = (1,1,...,1) € R* ofall ones. The latter modification 
makes sure that 
1. If the system 
H(w)w =v 


is solved for Aw € RX with v € Wx, then Aw € Wx automatically. 

2. If the eigenvalues of H(w), considered as an operator on Wx, are con- 
tained in [\_, à+], the eigenvalues of H(w) are contained in [\_, A+] 
as well. Furthermore, if |D;| > 0 for all i = 1,2,..., K, the function g 
is even strictly convex, and, thus, the eigenvalues of H(w) on Wx are 
strictly negative, i.e., A; < 0. As a consequence, also the eigenvalues 
of H(w) are strictly negative, ensuring invertibility of H(w). 
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These two properties are useful when implementing Newton’s method, 
see section 3.3.4. Item 1 means that we may essentially forget about the 
space Wx, and work in RX without second thoughts. The second item 
makes sure that our modification does not lead to an ill-conditioning of 
the linear system determining (3.24) the Newton increment. Indeed, for 
any A > 0, the matrix 

HALL” 


is negative semi-definite. If A is not chosen carefully, it may deteriorate 
the conditioning of the matrix in question. We know that 


tr(H(w)) 
— eE [à 

K € [ ’ +l, 
as it corresponds to the mean of the eigenvalues. Furthermore, the unit 
vector 1/vV K 1 spans the kernel of H(w) (at least on the set where g 
is strictly concave), and 1/K 11” is the orthogonal projector onto this 


kernel. Combining these observations leads to item 2. 


3.2.2 Generating centroidal Laguerre tessellations 


For a tessellation {D;} of Y denote by c; € Y a centroid of D;, i.e., a 
minimizer of 


ceY 


/ d(c, x)? de — min. 
D; 


The periodic boundary conditions may give rise to special situations if 
the cells D; "wrap around Y". As a simple example, for the single cell Y, 
any c € Y will be a centroid. 

However, if the Laguerre cell D; is sufficiently small, i.e., it may be 
displaced in Y to a configuration not touching the boundary OY CR, 
the classical representation as the center of mass 


sad 
Gi = — xdr 3.10 
Di Jp, on 
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(relative to a suitable local coordinate system) may be used. We shall 
rely upon the latter expression for practical purposes. 

A Laguerre tessellation is called centroidal (Brieden and Gritzmann, 
2012) if each seed is a centroid of its corresponding cell. Although, 
when looking at solidification of polycrystalline materials, there does 
not seem to be any requirement that the nucleation site corresponds to 
the grain centroid, it has been observed, see Bourne et al. (2020), that 
centroidal Laguerre tessellations exhibit a better shape regularity than 
general Laguerre tessellations. 

Suppose a set of (mutually different) seeds (x1, £2,..., £K) € YK and 
volume fractions ¢ € R* with ¢; > 0 and Sa i = 1 are given. By 
the reasoning of the previous section, we find a set of weights w € 
RX whose associated Laguerre tessellation with seeds (x1, £2,..., £x) 
has cells with volume fractions ¢. Denote by (ci,c2,...,ck) € Y* 
the centroids of the respective Laguerre cells. A simple method for 
iteratively approaching a centroidal Laguerre tessellation with volume 
fractions ¢ as above is to use the centroids as new seeds. A pseudo-code 
of this classical Lloyd-type algorithm (Bourne and Roper, 2015) is shown 
in Alg. 1. 


Algorithm 1 Lloyd-type centering (¢, maxit, tol) 


1: Determine initial guess (x?,79,..., 7%) of seeds 

2 k+-0 

3: residual + +00 

4: while k < maxit and residual > tol do 

5: compute Laguerre tessellation with seeds (x7, r§,..., xt) 
for volume fractions ¢ 

6: compute centroids (ch, ck,..., c% ) by (3.10) 

7: ee ee cea VS (Crane) 

8: ke k+1 

9: update residual, see equation (3.11) 


10: end while 
11: return £, residual 
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Several remarks are in order: 


1. The initial guess of seeds may be determined by pseudo- or quasi- 
random number generators. 


2. For computing the weights of the Laguerre tessellation for given 
seeds and target volume-fractions, any of the algorithms described in 
section 3.3 may be used. 


3. For assessing convergence at the k-th iteration, we use the residual 


2 
residual = = 1 Pda} c$) ‘ (3.11) 


max L;)? 


i.e., the L?-norm of the difference between seeds and centroids, 
weighted by volume-fraction and normalized by the maximum edge 
length of the unit cell. A typical value for tol is 1071. 

4. The Lloyd-type algorithm presented can be shown to converge by 
using a monotonicity argument (Bourne and Roper, 2015). 


The problem of obtaining a centroidal Laguerre tessellation may also 
be written as an energy-minimization problem, see Liu et al. (2016). 
However, the resulting optimization problem admits a multitude of 
minimizers, i.e., we face a non-convex optimization problem, in general, 
and it is difficult to design powerful computational strategies. For 
instance, there are examples where the BFGS method does not converge, 
even with exact line search, see Dai (2012). 

As an alternative, we investigate Anderson acceleration (Anderson, 
1965), which is a technique for accelerating general fixed point iterations 
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where, in our context, x” € Y* and the mapping 
F:Y5 >Y" (3.12) 


is given by computing the K centroids of the Laguerre tessellation, 
see Alg. 1. 

More precisely, the Anderson-accelerated version of the fixed point 
iteration (3.12) with depth m > 0 stores the last m iterates 


ar gel o. „hm 


and their images w.r.t. F 
Fe Fe Rz? 


computes a" € R™ as the minimizer of the problem 


2 


>, Gal Fie) — gtt) — min (3.13) 
i=1 = 
and proposes the next iterate as 
ght) = ara, (3.14) 
i=1 


Several remarks are in order. 


1. Problem (3.13) is a quadratic problem with linear constraint. In 
particular, its KKT conditions are described by an (m+ 1) x (m + 1) 
linear system. As m is typically small (not exceeding 10), this problem 
can solved directly with negligible computational overhead. 

2. Recall that we are using periodic boundary conditions for the cell 
Y = [0, L1) x [0,L2) x ... x [0, La). Thus, for any y € Y and z € Z, 
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there are countably many vectors € € R? with 
Yi = Zi + & mod Li (=1,2,...,d). (3.15) 


Indeed, for any such €, and any d-tuple (ki, ka,..., ka) of integers, we 
may form the vector 


d 
E=£+ y kiLiei, (3.16) 


i=l 


which also satisfies (3.15) with €; replaced by &. In turn, any two 
vectors € and £ satisfying equation (3.15) are related by (3.16) by a 
suitable tuple (kı,ka,...,ka) of integers. Thus, we define y — z as 
the vector € satisfying (3.15) with minimal length (which is well- 
defined if y and z are sufficiently close in the periodic distance). 
Using this definition K-times clarifies how the difference vectors 
F(a*+1—-*) — g*+1~* in (3.13) have to be understood. 

3. The system matrix for the KKT-system of problem (3.13) may be 
updated at each iteration by computing m inner products. 

4. For the first steps, i.e., where k < m, m = k — 1 is used for the current 
depth. We implement this by forcing the last few components of a” 
to zero. 

5. For linear equations, Anderson acceleration of depth m is "essentially 
equivalent" to GMRES(m) (Saad and Schultz, 1986), see Toth and 
Kelley (2015). 

6. For general nonlinear equations, Anderson acceleration may be in- 
terpreted as a Quasi-Newton method, see Fang and Saad (2009). 
More precisely, it builds up an approximate Hessian from the last m 
secants. Under a suitable non-degeneracy condition it can be shown 
rigorously that Anderson acceleration improves the convergence 
rate of linearly converging fixed point iterations (Evans et al., 2020). 


3.3 Optimization methods 


For degenerate problems, safeguarding strategies are necessary for 
improving convergence (Fu et al., 2020). 
7. The Anderson-accelerated Lloyd-type centering generalizes the re- 
laxed (or damped) Lloyd method for m = 1 with variable relaxation. 
The resulting algorithm is summarized in Alg. 2. For m = 0, we 
(formally) recover the original algorithm Alg. 1. For our numerical 
experiments, we set m = 4. 


Algorithm 2 Anderson-accelerated Lloyd-type centering 
(ġ, m, maxit, tol) 


: Determine initial guess (x), x9,...,x%) of seeds 
k-0 

: residual + +00 

: while k < maxit and residual > tol do 
compute Laguerre tessellation with seeds (x7, 2%, ..., 2%) 
for volume fractions ¢ 

compute centroids (cf, ch, ..., c#-) by (3.10) 
update the linear system for solving (3.13) 
determine the vector a” by solving (3.13) 
update at! + 0", aF (etti) 

10: k&k+1 

11: update residual, see equation (3.11) 

12: end while 

13: return £, residual 


akon H 


3.3 Optimization methods 


In this section, we describe pertinent optimization algorithms used for 
solving the concave optimization problem (3.4) 


— 
Kw) > ee 
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for the unknown weight vector, and given seeds and volume fractions 
(which we assume to be fixed throughout this entire section). To bring 
it into a more convenient form, we shall equivalently rewrite the latter 


problem as a minimization problem via 


i 3.17 
u) min. (3.17) 
with f(w) = —g(w), i.e., in terms of the convex function f. In particular, 


we shall use descent methods instead of ascent algorithms, for instance. 


3.3.1 The gradient descent method 


Recall that the gradient of the function f (3.17) is given, see equation 
(3.7), component-wise by 


Thus, a simple gradient descent method 
wet? = w” —a,Vi(w"), (3.19) 


involving a positive sequence of step sizes a,,, may be used for solving 
(3.17), as pioneered by Aurenhammer et al. (1998). If the (local) Lipschitz 
constant of the gradient of f is known to be M, the iterative scheme 
(3.19) can be shown to converge for the choice a,, = 1/M, see Nesterov 
(2004). However, for convex f, the convergence rate is only logarithmic, 
in general. If, in addition, f is m-strongly convex in the vicinity of the 
minimizer, the convergence rate improves to be linear, and the choice 
Qn = 2/(M + m) optimizes the convergence rate (estimate). Notice that 
the required iteration count is proportional to M/m. 

Unfortunately, for the problem (3.17) both the Lipschitz constant and the 
strong convexity constant are not easily available. Of course, both may be 
obtained as the minimum and maximum eigenvalue of the Hessian (3.8). 
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However, even computing the entries of the Hessian matrix (not even 
accounting for extracting the extremal eigenvalues) leads to significant 
computational overhead. 

We are only interested in the qualitative performance of the gradient 
scheme (3.19) in order to comparing it asymptotically to more elaborate 
minimization schemes. Thus, taking into account the necessary dimen- 
sion of the step-size (the weights have dimensions length”, the gradient 
has dimension d), we simply use 


an =0.1 (LiLo Lajt. (3.20) 


Aurenhammer et al. (1998) propose using a step-size of Polyak-type 
(Polyak, 1987). However, the latter requires a good estimate of the opti- 
mal objective-value (which is not available, in general) to be competitive. 
Thus, we stick to the simple choice (3.20). 


3.3.2 The Malitsky-Mishchenko solver 


Malitsky and Mishchenko (2019) designed an adaptive strategy for 
choosing the step size in gradient methods which is globally convergent 
for general convex continuously differentiable functions (not necessarily 
with globally Lipschitz-continuous gradient) that automatically exploits 
strong convexity. The method is non-monotone, but significantly outper- 
forms gradient descent methods with fixed step-size, in general. 

The Malitsky-Mishchenko method involves, in addition to the step-size 
an, another sequence {6,,} which is (formally) initialized by 0) = +x. 
Then, the scheme proceeds by updating a running estimate of the inverse 
Lipschitz constant of the gradient as the step-size 


_ lw” = wr 
an = min (VIF Tresen gan vA) 2D 
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and gradient descent steps (3.19), see Alg. 1 in Malitsky and Mishchenko 
(2019). The sequence 9, is updated in each step using the step-size 
On, = an /A,-1. The method requires an initial step size ag, which we 
take as for the gradient descent scheme (3.20). 


3.3.3 The Barzilai-Borwein method 


Originally devised as a special type of Quasi-Newton method, the 
scheme proposed by Barzilai and Borwein (1988) may also be interpreted 
as an adaptive gradient descent method using the step size 


|w” _ w ||? 


(Vf) = Var YP (wr — wy 


Qn = (3.22) 
For the Barzilai-Borwein method also take the step-size of the gradient 
descent scheme (3.20) as ao. 

The Barzilai-Borwein (3.22) leads to r-linear convergence when ap- 
plied to strongly convex quadratic objective functions (Dai and Liao, 
2002). By a perturbation argument (Liu and Dai, 2001), the method 
can also be shown to be locally q-linearly convergent when applied to 
general strongly convex and twice differentiable objective functions. 
However, there are counter-examples to global convergence for the 
latter class, leading to the development of non-monotone backtracking 
techniques (Fletcher, 2005) and step-size limiting approaches (Burdakov 
et al., 2019), for instance. When applying the Barzilai-Borwein method 
to computational micromechanics (Schneider, 2019), these globalization 
techniques were not necessary for the computational examples consid- 
ered. On the contrary, in most cases they decreased the performance 
profile. For computing Laguerre tessellations, we made a similar observa- 
tion, i.e., we did not encounter a single case where the Barzilai-Borwein 
method diverged, see section 3.4 below. Thus, we shall rely upon the 
Barzilai-Borwein prescription (3.22) without globalization, although we 
are not aware of a mathematical justification for this choice. 
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3.3.4 Newton’s method 


For finding a critical point of a twice continuously differentiable function 
f, Newton’s method iteratively updates 


wrth = w” ae Ant”, (3.23) 
where the Newton increment €” solves the equation 
V? rw") E = —VF(w"). (3.24) 


ay € (0, 1] is a damping factor (similar to the step-size of the gradient 
method). If w* is a critical point of f with non-degenerate Hessian 
V? f(w*) and V?f satisfies a Lipschitz condition close to w*, Newton’s 
method with a, = 1 can be shown to be quadratically convergent in 
the vicinity of w* (Kantorovich, 1948). However, global convergence is 
not achievable that way. Among the various globalization techniques 
for Newton methods (Nocedal and Wright, 1999), the damped Newton 
method (3.23) is particular popular. Far away from a critical point, 
Qn < 1 may be necessary to ensure global convergence. However, 
in the vicinity of a critical point, no damping is necessary (a, = 1 
for n > no). In this way, the local convergence speed of Newton's 
method is preserved. 

If f is (strictly) convex, €” determined by (3.24) is always a descent 
direction for f at w”. Thus, a proper backtracking strategy allows us to 
determine a suitable step size a,,. More precisely, given a backtracking 
factor p € (0, 1), starting from aù, = 1, the step size is iteratively shrunk, 


k+1 


k+1 — pak, until actual descent is observed, either in function value 


i.e., a 


Fw” + ane”) < f(w") 
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or in its gradient 
VW" + one" )IL < IV E(w"). 


For the problem at hand (3.17), we shall rely upon the latter criterion, 
as V f is cheaper to compute than f. Under suitable conditions, this 
backtracking strategy can be shown to be globally convergent under 
technical conditions, see Boyd and Vandenberghe (2004). We take 
p = 1/2 as our backtracking factor. 

In the context of Laguerre tessellations with prescribed volumes, a 
further issue needs to be addressed. For computing the Newton in- 
crement €", the equation (3.24) needs to be solvable. In section 3.2.1, we 
pointed out that constant vectors are in the kernel of V? f(w) for any 
w € R. However, this is no problem if equation (3.24) is interpreted as 
an equation on Wx, i.e., Vf(w") € Wk is given and ¿” € Wx is sought. 
A simple remedy consists in adding the (suitably scaled) projector onto 
the kernel onto V? f(w”), see equation (3.9). 

Even restricted to Wx, the function f is not globally strictly convex. At 
such a point w, also the Hessian will be degenerate (even when restricted 
to Wg). More precisely, taking a look at the Hessian of g, we see that 
an empty Laguerre cell D; implies that both the i-th row and the i-th 
column only contain zeros. 

A workaround has been proposed by Lévy (2015) by adding the con- 
straint of non-empty cells to the backtracking strategy of Newton’s 
scheme. This approach has been analyzed in Kitagawa et al. (2019) and 
shown to be globally convergent when started at a feasible point (i.e., 
with non-empty Laguerre cells). 

However, for grain-size distributions with large variation, we found 
this cautious backtracking strategy to be computationally inefficient. 
Although the method is quadratically convergent locally, the first few 
steps may require extensive backtracking. As a simple alternative, we 
add - in case of vanishing i-th cell - a suitably scaled 1 on the i-th diagonal 
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entry of the Hessian. In this way, we ensure obtaining a descent direction 
and avoid the extensive backtracking. 


3.3.5 BFGS 


For large number K of grains, the computationally most demanding 
step of Newton’s method (3.23) is solving the linear system (3.24) which 
determines the Newton increment. If a direct solver is used, on the order 
of K? floating point operations are required, which may be excessive for 
large K. 

To reduce the effort, Quasi- Newton methods approximate the inverse 
of the Hessian V f(x)! directly, see Nocedal and Wright (1999). The 
most popular Quasi-Newton method is the Broyden-Fletcher-Goldfarb- 
Shanno method (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 
1970), or BFGS, in short. In a nutshell, for every n, the next iterate is 
computed as for Newton’s method (3.23) 


wrth = w” Az Ané” 
but the increment is determined by 
er = -HV fu"), 


where H, € R*** is a sequence of approximations to the inverse 
Hessian, which is, for BFGS, updated via 


Cae Ir ye Hn Yn T 1 
(sin)? ee EY 


Ana, = Hn + (Anns, + Sny Hn) 


(3.25) 


where yn = Vf(w"t!) — Vf(w") and sn = w”+! — w” . For strictly 
convex objective function f, the BFGS update (3.25) preserves symmetry 
and positive definiteness of the Hessian approximation, see Nocedal and 
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Wright (1999). Thus, if the initial approximation Ho is symmetric positive 
definite, all further H,, will remain symmetric and positive definite. 
The effectiveness of BFGS crucially depends on the initial choice Hp. A 
common choice (Nocedal and Wright, 1999) is to use the Barzilai-Borwein 
step size (3.22), i.e., 


E lw? — w~]? 
Ho = pa) = Vf) A 


provided some inital gradient step has been taken, as for the Barzilai- 
Borwein method. When enriched by a backtracking strategy for 
determining the step size an in (3.23), BFGS can be shown to converge 
for strongly convex functions (and further technical conditions) with a 
superlinear rate (Nocedal and Wright, 1999). 

Notice that we do not need to pay specific attention to the kernel of 
V? f consisting of component-wise equal weights. This is a result of H„ 
approximating the inverse of V?f on Wr. Indeed, Ho is a multiple of 
the identity (with a proper scaling) and the BFGS update (3.25) only 
modifies Wg. 

The BFGS method reduces the complexity required for the New- 
ton step from O(K?) to O(K?) (due to the update of H,, and the 
matrix-vector product required for computing €"). To further re- 
duce the complexity to O(K), limited-memory BFGS (L — BFGS), 
see Nocedal (1980), may be used. In general, however, the super- 
linear convergence is lost. Also, L — BFGS is often inferior to the 
Barzilai-Borwein method, as the latter scheme’s non-monotonicity 
can be advantageous (Schneider, 2019; Wicht et al., 2020b). 
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3.4 Numerical demonstrations 


3.4.1 Setup and used hardware 


The algorithms described in sections 3.2.2 and 3.3 were implemented in 
Python 3.7 with Cython extensions. For computing Laguerre tessella- 
tions, we rely upon the Voro++ (Rycroft, 2009) code which we exposed 
to Python via Cython bindings using an OpenMP parallelization. 

To reduce the stochastic influence, we use the Sobol low-discrepancy 
sequence (Sobol, 1967) for drawing log-normal grain-size distribu- 
tions (Johnson et al., 1994). The numerical experiments were conducted 
on desktop computer with 6 3.7 GHz CPUs and 32 GB RAM. 

For solving the linear system for the Newton increment (3.24), we rely 
upon NumPy’s linalg.solve routine, which wraps LAPACK (Anderson 
et al., 1999). The "higher-level" vector and matrix operations are realized 
by NumPy, for instance the BFGS update (3.25). 

For all convex optimization methods of section 3.3, we initialize the 
weights to zero. Using the estimate of Lyckegaard et al. Lyckegaard 
et al. (2011), as suggested in Bourne et al. Bourne et al. (2020), did not 
decrease the run-time, on average. 

To ensure comparability between different solvers, all experiments are 
initialized with identical seeds. We realize this by using the 3D Sobol 
sequence on the (scaled) unit cube. The first few elements of the 3D Sobol 
sequence are a little pathological, as they correspond to the corners of the 
cell and some other very regularly positioned spots, leading to extremely 
regular Laguerre tessellations. Therefore, we skip a few elements of the 
Sobol sequence. To ensure comparability, we fix the number of elements 
to be skipped as 1234 (for no specific reason). 

At this point, we wish to draw attention to the convergence criteria 
and tolerances used. For the convex optimization solvers, we test 
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convergence by 
voll <tol (3.26) 
minj~ di 
with tol = 10~?. This choice ensures that the volume fraction of the 
smallest grain is accurate to two significant digits (which we assume 
sufficient for materials-science purposes). For centering, we use equation 


(3.11), i.e., 


Ko, dak. ck)2 
= i d(x}, c7) 2 tol 


(max L;)? = 


with tol is 107+. As Bourne et al. (2020), we noticed no visible improve- 
ment of the tessellation quality for lower values of tol. In contrast, the 
increase in computation times became quite visible. 

We abbreviate grain-size distribution as GSD throughout section 3.4. For 
all the experiments in this section, we use cubical cells, i.e., Li = La = Ls. 
The prescribed algorithms work for anisotropic unit cells ("bricks") 
as well. 


3.4.2 Computing the weights of a Laguerre tessellation 
with prescribed volume fractions 


The overall goal of this section is to evaluate the solvers introduced in 
section 3.3. For that purpose, we evaluate their iteration counts and 
run-times for problems of different complexity and grain count. 

As our first example, we seek Laguerre tessellations with a mono-sized 
GSD, i.e., for K grains, each grain should have a volume fraction 1/K. 
As our smallest number of grains, we choose K = 250. For smaller 
grain count, the run-times become even more negligible and a sensible 
evaluation is difficult. Starting from K = 250, we subsequently double 
the grain count up to 16000 cells. 
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Larger numbers are possible (and computationally feasible), but prob- 
ably limited in terms of applicability to computational homogeniza- 
tion (Shantraj et al., 2015; Wicht et al., 2020a). 

The iteration counts of the investigated solvers are listed in Tab. 3.1. The 
gradient method (3.19) serves as our basic reference. Even for K = 250, 
more than 300 iterations are required. For higher K, the iteration counts 
are proportional to K, i.e., they roughly double when doubling the 
number of grains. 

For the Malitsky-Mishchenko method (3.21), the iteration counts are 
significantly lower than for the gradient method, roughly by one or 
two orders of magnitude. Still, we see a dependence of the iteration 
count on the grain number. More precisely, the iteration count depends 
approximately linearly on K, but with a non-trivial offset at K = 0. 
The Barzilai-Borwein method leads to a consistently lower iteration 
count than the Malitsky-Mishchenko method, except for K = 250, but 
no more than a factor of 2. BFGS (3.25) outperforms the Barzilai-Borwein 
scheme consistently in terms of iteration counts, also by no more than a 
factor of two. 

Newton’s method, see section 3.3.4, operates on a completely different 
level in terms of iteration counts, requiring 2 — 3 iterations for all grain 
counts considered. In particular, we see quadratic convergence in action. 


number of grains K 250 500 1000 2000 4000 8000 16000 


Gradient method 352 628 1017 1817 3275 5411 9758 
Malitsky-Mishchenko 16 28 37 40 62 87 124 
Barzilai-Borwein 18 20 31 32 39 59 66 
BFGS 11 15 18 19 22 32 38 
Newton 2 3 3 3 3 3 3 


Table 3.1: Iteration counts for different solvers for mono-sized GSD. 
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When evaluating computational methods, apart from the iteration 
counts, also the run-times are of particular importance, see Tab. 3.2. 
Even at first glance we see that the gradient method is not competitive 
in terms of run-time, which is certainly not surprising. The other four 
solvers perform quite well in terms of run-time for the majority of 
problems, requiring only a few seconds. 

The Malitsky-Mishchenko and the Barzilai-Borwein method are both 
gradient methods of similar type, i.e., the computational expense of their 
individual steps are comparative (see also Tab. 3.3). For K < 2000, both 
of their run-times are similar (and similarly negligible). Only starting at 
K = 4000, differences are visible. For K = 16000, the Barzilai-Borwein 
method is about twice faster than Malitsky-Mishchenko. 

Jumping ahead to Newton’s scheme for the time being, we see that 
for small K (up to K = 2000), Newton’s method is the fastest method 
investigated. For larger number of grains, Newton’s scheme falls back 
behind the Barzilai-Borwein method in terms of run-time. 

BFGS appears to be in between Barzilai-Borwein and Newton. For small 
grain count, it is slower than Newton’s method, and for large K, it is 
slower than the Barzilai-Borwein scheme. 


number of grains K 250 500 1000 2000 4000 8000 16000 


Gradient method 1.59 3.98 9.63 30.54 100.73 325.34 1155.00 
Malitsky-Mishchenko 0.18 0.28 0.47 0.77 2.02 5.35 14.69 
Barzilai-Borwein 0.19 0.27 0.54 0.65 1.30 3.69 8.21 
BFGS 0.16 0.27 0.50 1.29 5.31 26.93 103.64 
Newton 0.14 0.16 0.20 045 1.52 5.41 32.67 


Table 3.2: run-time in s for different solvers for mono-sized GSD. 


More insight into these observations may be gained by investigating 
the run-time per iteration of the different solvers under consideration, 
see Tab. 3.3. Remarkably, the three gradient-type methods require 
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roughly the same computational effort per iteration. The expense is 
slightly higher for the Barzilai-Borwein method than for the other two 
gradient-type methods because the Python overhead plays a more 
distinct role, essentially due to the lower number of iterations required 
for Barzilai-Borwein method. Recall that we exposed the C++ - library 
Voro++ (Rycroft, 2009) to Python via Cython, and the major computa- 
tional workload of the gradient-type methods amounts to computing 
Laguerre tessellations via Voro++. This library computes the Laguerre 
cells individually (and sequentially), using cell-linked lists (Allen and 
Tildesley, 1987). As our initial seeds are more or less equally distributed, 
and - for mono-sized GSD - we look for cells of equal volume, the effort 
for computing a single Laguerre cell is effectively independent of the cell 
count. In turn, the total computational expense for computing the entire 
Laguerre tessellation is proportional to K. We see this trend quite clearly 
in Tab. 3.3 for the gradient-type solvers, highlighting the quality of the 
Voro++ - implementation. Notice the non-perfect scaling of the run-times, 
which is a result of the Python overhead and becomes less pronounced 
for large K. A close look at the run-time per iteration for BFGS reveals 
a quadratic dependence on K. In contrast, the effort required for a 
Newton step grows cubically in K. These observations also clarify why 
the overall run-times of BFGS lag behind those of Newton’s method for 
large K. BFGS has a lower effort per iteration than Newton's method, but 
the overall iteration count is also significantly larger than for Newton’s 
method, tipping the overall scale towards Newton. 


number of grains K 250 500 1000 2000 4000 8000 16000 


Gradient method 0.45 0.63 0.95 1.68 3.08 6.01 11.84 
Malitsky-Mishchenko 1.12 0.99 1.27 1.93 3.27 6.15 11.84 
Barzilai-Borwein 1.06 1.35 1.75 2.02 3.34 6.26 12.45 
BFGS 1.48 1.79 2.78 6.79 24.15 84.14 272.73 
Newton 6.92 5.17 6.80 14.86 50.58 180.18 1088.62 


Table 3.3: run-time per iteration in 10”2s for different solvers and a mono-sized GSD. 
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To finish this first example, we closer examined the convergence behavior 
of the solvers under consideration. For that purpose, for K = 1000 
grains, we solve equation (3.2) up to a precision of tol = 10~'°. This 
level of accuracy is certainly beyond what is necessary for applications, 
but gives interesting insights into the internal mechanisms of the solvers. 
The residual (3.26), is plotted vs. the iteration count in Fig. 3.2. We 
did not include Newton’s method due to its quadratic convergence. It 
appears that all shown solvers converge linearly. The gradient method 
converges slowly, but in a steady and monotonic fashion. Both the 
Malitsky-Mishchenko (MM) and the Barzilai-Borwein (BB) method do 
not exhibit monotone convergence behavior but lead to an oscillation 
of the residual by one or two orders of magnitude. Overall, both 
schemes converge r-linearly. This complies with theoretical estimates, 
see section 3.3. Overall, the convergence speed of BB is higher than 
for MM. BFGS also converges linearly, but in a monotonic fashion. 
Put differently, we have q-linear convergence. BFGS exhibits a higher 
convergence rate than the other solvers shown. However, no superlinear 
convergence is observed. This is consistent to theory, as to enable 
superlinear convergence, the approximation to the inverse Hessian needs 
to be close, as well, see section 8.4 in Nocedal and Wright (1999). As 
the Hessian is a 1000 x 1000-matrix, the less than 100 iterations (and in 
turn, the corresponding matrix updates (3.25)) are simply insufficient to 
match the true inverse Hessian closely. 
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Figure 3.2: Convergence behavior of different inner solvers for the mono-sized GSD and 
K = 250 grains. 


As our next computational setup we consider a grain-size distribution 
with exactly two (distinct) volume fractions, each of which is equally 
probable. This may be viewed as a simplified bimodal grain-size distribu- 
tion (Flipon et al., 2020a). Following Bourne et al. (2020), we use a factor 
of 5 between the two respective volume fractions. Thus, we consider K 
grains with volume fraction 1/3K and K grains with volume fraction 
5/3K, each. To keep consistency to the mono-sized GSD, the total 
number of grains 2K is varied between 250 and 16000, incrementally by 
factors of 2. Due to its lack of competitiveness for the mono-sized GSD, 
we did not consider the gradient descent scheme any further. 


number of grains K per phase 125 250 500 1000 2000 4000 8000 


Malitsky-Mishchenko 34 41 86 83 110 203 295 
Barzilai-Borwein 25 32 40 48 63 71 114 
BFGS 19 20 29 28 32 41 59 
Newton 3 5 6 4 4 5 7 


Table 3.4: Iteration counts for different solvers and the two-sized GSD. 
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The iteration counts for the four investigated solution methods are listed 
in Tab. 3.4. Both for Malitsky-Mishchenko and for Barzilai-Borwein, 
the iteration count increases monotonically with K. A similar trend 
may be observed for Newton’s method and BFGS. However, there are 
exceptions between K = 250 and K = 1000. It is also evident that the 
iterations counts for the two-sized GSD are consistently higher than for 
the mono-sized GSD, see Tab. 3.1. 

Comparing the iteration counts of Malitsky-Mishchenko and Barzilai- 
Borwein, we see that for increasing K, the discrepancy between the 
two solvers, in terms of iteration count, is increasing. For K = 
8000, Malitsky-Mishchenko requires about three times the iteration 
count of Barzilai-Borwein. 

As for the mono-sized GSD, BFGS requires less iterations than Barzilai- 
Borwein, roughly by a factor of two for the larger K. 

For Newton’s method, no quadratic convergence behavior can be 
observed for the two-sized GSD. Indeed, up to 7 Newton iterations are 
necessary to trigger the convergence criterion. Intuitively, this appears 
clear. Newton’s method exhibits its strength close to the minimum. 
However, the seeds are evenly distributed, and the weights are equal 
initially. Thus, initially, the computed Laguerre cells will have roughly 
equal volume. In particular, they do not appear to serve as a promising 
starting value for Newton’s method. Therefore, a few preliminary 
descent steps are necessary to bring the iterates into the domain of 
quadratic convergence. 


number of grains K per phase 125 250 500 1000 2000 4000 8000 


Malitsky-Mishchenko 0.25 0.37 0.97 152 3.51 12.73 36.55 
Barzilai-Borwein 0.24 0.32 053 0.93 2.08 4.44 14.51 
BFGS 0.22 0.29 0.90 1.96 7.86 34.08 171.39 
Newton 0.14 0.19 0.27 0.52 1.85 871 92.74 


Table 3.5: run-time in s for different solvers and the two-sized GSD. 
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The required run-times are listed in Tab. 3.5. In general, the run-times 
exceed those of the mono-sized problem, see Tab. 3.2. The Malitsky- 
Mishchenko method consistently lags behind the Barzilai-Borwein solver 
in terms of run-time. Newton’s method is very strong except for the 
highest grain size considered. More precisely, up to K = 2000, Newton’s 
method outperforms all other solvers considered. For K = 4000 and 
K = 8000, the Barzilai-Borwein method is fastest. Again for the two- 
sized scenario, BFGS does not perform admirably. 
Notice that up to K = 2000, the run-times of Barzilai-Borwein and 
Newton’s method are minimal (up to about 2s). Also, for the highest 
grain count considered, 16000 grains, the Barzilai-Borwein solver needs 
less than 15s. Compared to the mono-sized problem, see Tab. 3.2, this 
amounts to an overall increase by a factor less than two. 
For further analysis, we discard BFGS. 
To model more realistic grain-size distributions, we use a log-normal 
distribution for the equivalent radius of the grains. To be more precise, 
for a grain of volume V, we associate the radius r of a sphere with 
precisely this volume to it 
V= a5 i \ Ln 
37° 
It is empirically known (Spettl et al., 2014) that the distribution of this 
equivalent radius follows a log-normal distribution with probability 
density function (Johnson et al., 1994) 


p(r)=-- = exp ( gr = n) ) ; (3.27) 


depending on the two real parameters u and o. The mean and standard 
deviation compute as (Johnson et al., 1994) 


12 1,2 
mean = e"+27 and stdev = ett3” Ver? — 1. 
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As the mean equivalent radius plays no role in the geometry generation, 
we set it equal to unity. In turn, we obtain the relation u = -1/20°, 
reducing the original two parameters to the single parameter ø. In turn, 
the expression for the standard deviation simplifies to 


stdev = Ve?” — 1. (3.28) 


Thus, for given standard deviation stdev, we may determine the 
parameter o characterizing the log-normal distribution (3.27) with 


co = \/log(1 + stdev”). 


Thus, we may parameterize grain-size distributions by the standard 


mean | explicitly 


deviation. For vanishing standard deviation, we recover the mono-sized 
distribution. Spettl et al. (2014) report stdev = 0.35 as typical for a 
polycrystalline material that underwent grain-boundary migration by 
capillary effects. 
To sum up, we prescribe the unit cell size (L4, La, L3), the number K of 
grains and the standard deviation stdev of the log-normal distribution. 
Then, we sample radii 71, ra,...,rx according to the log-normal distri- 
bution by use of the Sobol sequence (Sobol, 1967). Then, we compute 
the appropriate volume fractions by 
= 
pie 1 r3 ' 
which serves as input for our optimization routines. 
For the experiment at hand, we fix the grain count to K = 1000 and 
vary the standard deviation from 0 to 0.5 in 0.05-steps. Fig. 3.3 contains 
images for microstructures generated for varying standard deviation. 
The seeds are identical for all structures. We see that some grains 
grow, whereas others shrink, for increasing standard deviation. Please 
note that the colors appearing in the figures of this article showing 
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polycrystalline microstructures are chosen randomly. Their sole purpose 
is to distinguish the individual cells. In particular, they are devoid of 
any physical meaning. 


(a) (c) 


(d) (e) (£) 


Figure 3.3: Influence of increasing standard deviation for polycrystalline microstructures 
with log-normally distributed equivalent radii and 1000 grains. The seed locations are 
fixed and the plots show different draws from the log-normal distribution. (a) stdev = 0.00, 
(b) stdev = 0.10, (c) stdev = 0.20, (d) stdev = 0.30, (e) stdev = 0.40, (f) stdev = 0.50. 


The iteration counts and run-times for increasing standard deviation 
and the Malitsky-Mishchenko and the Barzilai-Borwein scheme as well 
as Newton’s method are listed in Tab. 3.6. 
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MM BB Newton MM BB Newton 
0.00 37 31 3 0.35 0.30 0.11 
0.05 94 60 3 0.99 0.60 0.12 
0.10 134 69 3 1.41 0.75 0.13 
0.15 162 99 4 1.83 1.17 0.14 

5 
9 
5 


0.20 215 115 2.61 1.39 0.32 
0.25 272 137 3.48 1.74 0.58 
0.30 343 176 4.62 2.37 0.19 


0.35 460 168 5 7.03 2.42 0.19 
0.40 628 250 11 9.74 3.95 0.75 
0.45 934 359 25 15.07 5.96 2.93 
0.50 1742 499 48 29.77 8.85 6.05 


Table 3.6: Iteration counts and run-times for different solvers and varying standard 
deviation for K = 1000 cells. 


Apparently, as a general rule of thumb, increasing standard deviation 
also leads to an increase in iteration count and run-time for all solvers 
considered. This is not only caused by the increased difficulty, but 
also by the convergence criterion (3.26) in use. Indeed, the smallest 
among the volume fractions enters the convergence criterion. As the 
smallest volume fraction decreases for increasing standard deviation, 
the convergence criterion becomes stricter, as well. At this point, we may 
question the convergence criterion we use. However, this convergence 
criterion is standard, see, for instance, Lévy (2015). 

Taking a closer look at the iteration counts in Tab. 3.6, we see that the 
Malitsky-Mishchenko method consistently requires more iterations than 
the Barzilai-Borwein method. For stdev = 0.5, it even requires more than 
thrice the iteration count of Barzilai-Borwein. The required iterations 
for the Barzilai-Borwein method increases superlinearly in the standard 
deviation. It appears that the iteration count roughly doubles every 
0.1-step in standard deviation. 

For increasing standard deviation, Newton’s method is still superior, but 
with a smaller margin, because it can profit less from its local quadratic 
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convergence behavior. For low standard deviation (below 0.15), three 
iterations suffice. For intermediate standard deviation (less than 0.4), 
only a few Newton steps need to be taken before entering the domain 
of quadratic convergence. For larger standard deviation, the iteration 
counts increase significantly, rising up to 48 for stdev = 0.5. 

Recall that we used K = 1000 grains. For stdev = 0, we recover the 
mono-sized experiment. The run-times for this case lie well below one 
second. For Malitsky-Mishchenko, compared to the mono-sized case, 
for stdev = 0.5 the run-time is increased by a factor of 100. For the 
other two solvers, the increase is less severe, but still noticeable. Overall, 
for standard deviations up to 0.4, Newton’s method is extremely fast, 
with run-times below one second. Only for higher standard deviation, 
a significant increase in run-time is observed. The run-times for the 
Barzilai-Borwein method increase more steadily with increasing stan- 
dard deviation, and are only slightly above those of Newton’s method. 
Recall, from the previous experiments, that the optimal solver strongly 
depends on the number of grains present in the tessellation. For K = 
1000 grains, Newton’s method always proved the fastest. For the last 
experiment, we saw that the standard deviation of the grain-size distri- 
bution also plays a significant role when evaluating the performance of 
the solvers in question. 


3.4.3 Computing centroidal Laguerre tessellations for 
prescribed volume fractions 


In the previous section 3.4.2, we computed periodic Laguerre tessella- 
tions for prescribed volume fractions and fixed seeds. In this section, 
we wish to determine centroidal Laguerre tessellations for prescribed 
volume fractions, and add the centering algorithms of section 3.2.2 on top 
of the optimization algorithms of section 3.3. We have seen that either 
Newton's method or the Barzilai Borwein scheme proved computation- 
ally most efficient for the problems in section 3.4.2. Thus, we restrict 
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to these two solvers for the subsequent investigations. For the entire 
section, we use the convergence criterion (3.11) with tolerance tol = 107+. 
Recall that we check convergence for the convex optimization solvers by 
condition (3.26) with tol = 107?. 

Consistent to section 3.4.2, we first investigate a mono-sized GSD with 
increasing number of grains. 

The number of centering iteration counts is listed in Tab. 3.7. We see 
that the number of centering iterations depends only on the "outer" 
algorithm used, and is independent on the "inner" algorithm. This 
fact appears clear at first, but becomes more interesting as we chose 
the tolerance for inner residual (3.26) comparatively large. Indeed, 
for inexact Newton methods (Dembo et al., 1982), the choice of the 
inner tolerance has a strong influence on the overall performance of the 
scheme, see Dembo et al. (1982). However, this seems not to be the case 
for the problem at hand. 


number of grains K 250 500 1000 2000 4000 8000 16000 


BB + Alg. 1 23 19 16 12 10 7 6 
Newton + Alg. 1 23 19 16 12 10 7 6 
BB + Alg. 2 13 11 10 8 7 6 5 
Newton + Alg. 2 13 11 10 8 7 6 5 


Table 3.7: Centering iteration counts for different optimization solvers, combined with 
the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the 
mono-sized GSD. 


Comparing the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated 
version Alg. 2, we see Anderson acceleration decreases the iteration 
count in any case. However, for high grain count, the improvement is 
not large. The associated run-times are listed in Tab. 3.8. Up to K = 2000 
grains, Newton’s method is faster. For higher grain count, the Barzilai- 
Borwein method leads to a lower run-time. Comparing the centering 
algorithms, Anderson acceleration consistently lowers the run-times. 
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number of grains K 250 500 1000 2000 4000 8000 16000 
BB + Alg. 1 0.59 0.81 1.33 2.25 449 9.61 19.27 
Newton + Alg. 1 0.32 0.42 073 185 7.51 23.81 140.84 
BB + Alg. 2 0.49 0.72 1.24 1.98 3.90 9.06 18.03 
Newton + Alg. 2 0.28 0.33 0.55 1.39 648 23.66 131.66 


Table 3.8: Run-time in s for different optimization solvers, combined with the Lloyd-type 
algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the mono-sized GSD. 


More insight might be gained from Fig. 3.4, where the number of inner 


iterations is shown vs. the current centering iteration for K = 250. For 
Lloyd’s method, see Alg. 1, the inner iteration count decreases monoton- 
ically for increasing centering iteration. This property is preserved by 


Anderson acceleration (except for a single centering iteration). 


To sum up, for large grain counts, i.e., K > 4000, combining the Barzilai- 
Borwein method with Alg. 2 performs best, requiring less than 20s for 


the largest problem considered. 
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— BB+Alg.1 — BB+Alg.2 
— Newton + Alg. 1 — Newton + Alg. 2 
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Figure 3.4: Inner iterations vs centering iterations for the outer solvers for the mono-sized 
GSD and K = 250 grains. 


As our next example, as in section 3.4.2, we consider a two-sized grain- 
size distribution with a ratio of 5 between the volume fractions of the 
two "phases", each of which has K grains. 

The centering iteration counts are listed in Tab. 3.9. Compared to the 
mono-sized GSD, the iteration counts are slightly higher, on average. 
Still, the observations for the mono-sized GSD also apply for this two- 
sized case. 


number of grains K per phase 125 250 500 1000 2000 4000 8000 


BB + Alg. 1 31 20 17 13 9 i 6 
Newton + Alg. 1 31 20 17 13 9 K 6 
BB + Alg. 2 18 13 12 10 7 6 5 
Newton + Alg. 2 18 13 12 10 7 6 5 


Table 3.9: Iteration counts for different optimization solvers, combined with the Lloyd-type 
algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the two-sized GSD. 
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As already observed in Bourne et al. (2020), the time per centering 
iteration decreases. Therefore, we investigate the corresponding run- 
times, contained in Tab. 3.10. Already for K = 4000, Barzilai-Borwein 
combined with Alg. 2 outperforms all other variants. For the highest 
grain count considered, the latter combination converges in about half a 
minute, which is about 33% higher than for the mono-sized GSD. 


number of grains K per phase 125 250 500 1000 2000 4000 8000 


BB + Alg. 1 140 1.57 2.42 404 642 12.45 33.41 
Newton + Alg. 1 0.51 0.63 1.34 3.05 16.01 55.08 406.16 
BB + Alg. 2 1.08 1.20 2.05 344 5.84 11.84 31.01 
Newton + Alg. 2 0.35 0.49 1.39 3.00 14.74 51.67 389.14 


Table 3.10: run-time in s for different optimization solvers, combined with the Lloyd-type 
algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the two-sized GSD. 


As our last example in this paragraph, we consider the problem of 
computing a centroidal periodic Laguerre tessellation with K = 1000 
grains with log-normally distributed GSD, cf equation (3.27), and vary- 
ing standard deviation. 

The number of centering iterations, see Tab. 3.11, is almost independent 
of the inner solver (different only by 1, at most). We see that the 
number of centering iterations increases monotonically with increasing 
standard deviation (except for Anderson and the transition from 0.45 
to 0.5). However, in contrast to the inner iterations, see Tab. 3.6, the 
increase is modest. For both centering algorithms, the run-times are 
comparable, and thus mainly determined by the inner solver. Taking 
a closer look at the run-times for Lloyd-type centering, we see that 
Newton’s method serves as the faster inner solver up to a standard 
deviation of 0.2. However, Barzilai-Borwein is only slightly faster. 
Indeed, for stdev = 0.5, both inner solvers lead to an overall run-time 
less than half a minute. 
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BB Newton BB Newton BB Newton BB Newton 

Alg.1 Alg.1 Alg.2 Alg.2 Alg.1 Alg.1 Alg.2 Alg.2 
0.00 16 6 10 10 1.37 0.59 1.06 0.45 
0.05 16 6 10 10 1.87 0.59 1.62 0.55 
0.10 16 6 10 10 2.35 0.69 1.99 0.48 
0.15 16 6 10 10 3.23 1.15 2.75 1.13 
0.20 16 6 10 10 4.09 3.35 3.41 3.47 
0.25 i T 11 11 4.92 8.07 4.37 8.36 
0.30 18 8 12 12 6.87 8.30 6.32 8.22 
0.35 18 8 12 13 8.68 9.67 7.63 10.12 
0.40 20 20 15 15 12.70 20.74 11.38 19.62 
0.45 21 21 17 17 18.53 25.66 17.61 27.22 
0.50 22 22 16 16 25.15 29.56 23.44 29.91 


Table 3.11: Centering iteration counts and run-times for different optimization solvers, 
combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version 
Alg. 2, and varying standard deviation. * means that the method did not converge in a 
reasonable amount of time. 


The generated centroidal polycrystalline microstructures, for different 
values of the standard deviation, are shown in Fig. 3.5. For comparison, 
the non-centroidal versions may be of interest, see Fig. 3.3. Recall that we 
used identical seeds for all cases, i.e., independent of standard deviation 
and whether centering is applied or not. Also, the coloring schemes 
are kept consistent. Thus, it is possible to match the various grains 
for all cases considered. Apparently, the centroidal tessellations in 
Fig. 3.5 exhibit a much higher shape regularity than the non-centered 
tessellations of Fig. 3.3. Also, the sphericity of the large grains becomes 
more pronounced for increasing standard deviation. 
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(d) (e) (£ 


Figure 3.5: Influence of increasing standard deviation for centered polycrystalline 
microstructures with log-normally distributed equivalent radii and 1000 grains. 
(a) stdev = 0.00, (b) stdev = 0.10, (c) stdev = 0.20, (d) stdev = 0.30, (e) stdev = 0.40, 
(£) stdev = 0.50. 
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3.5 Examples with a higher 
level of complexity 


After thoroughly assessing the performance of the convex-problem 
solvers and the centering, we shall investigate examples with a higher 
degree of complexity. 

We used the Barzilai-Borwein method to generate a microstructure with 
100000 = 10° grains of identical volume, up to an accuracy of 10~?, as 
before. For quasirandom initial positions, the Barzilai-Borwein method 
needed 282.85s, i.e., less than 5 minutes. The resulting microstructure 
is shown in Fig. 3.6a. Coupled to the Anderson-accelerated Lloyd 
algorithm, solved up to an accuracy of 1074, took 338.32s, i.e., less 
than six minutes. The resulting microstructure is shown in Fig. 3.6b 
and may be compared to its non-centered cousin, cf Fig. 3.6a. Even 
without the multi-level strategies of Mérigot (2011) and Lévy (2015), we 
are able to generate microstructures with high complexity to be used in 
crystal-plasticity FEM (Becker, 1991; Barbe et al., 2001; Roters et al., 2010) 
or FFT (Lebensohn and Rollett, 2020; Shantraj et al., 2015; Wicht et al., 
2020a) frameworks in a matter of seconds to a few minutes. 
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(a) (b) 


Figure 3.6: Generated microstructure with 100000 grains of identical volume. (a) Non- 
centroidal Laguerre tessellation, (b) Centroidal Laguerre tessellation. 


The framework presented may also be used for generating microstruc- 
tures with a morphological texture, i.e., a prescribed spatial gradient 
in the grain sizes. Such microstructures may be the result of high 
temperature gradients that may occur inside of polycrystalline metals, 
for example during welding (Lehto et al., 2016). 

We follow the ideas of Bourne et al. (2020) and restrict the initial positions 
of the seeds to particular subsets. Of course, the resulting grains need not 
be entirely contained in the pre-selected areas, in particular if coupled to 
the centering methods. Still, interesting microstructure emerge, as we 
shall demonstrate by an explicit example. 
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We consider a cubical cell [0, L] x [0,Z] x [0, L]. We shall divide the 
grains into three classes: 


1. 1024 grains with initial positions restricted to [0, 4] x [0, L] x [0, Z] 
accounting for 1/3 of the volume fraction. 

2. 256 grains with initial positions restricted to [£, 4] x [0, L] x [0, L] 
accounting for 1/3 of the volume fraction. 


3. 64 grains with initial positions restricted to [£,L] x [0,L] x [0,L] 


accounting for 1/3 of the volume fraction. 


Thus, we have three types of grains. Each type has the same total volume. 
However, the grain volume of class one is four times smaller than for the 
second class, and sixteen times smaller than the third class. Also, the first 
class is restricted to the first quarter of the x-face, whereas the second 
class is restricted to the second quarter, and the third phase makes up 
the last half of the x-axes. 

Generating the microstructure with Barzilai-Borwein and Anderson- 
Lloyd (with accuracy as above) took 2.0s. The resulting microstruc- 
ture is shown in Fig. 3.7b. The first class is shown in red-to-black, 
the second class in yellow, whereas the third class is characterized 
by a grayish coloring. 
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(a) (b) 


Figure 3.7: Generated polycrystalline microstructure with spatial size-gradient. (a) Slice 
through the center, (b) Volumetric view. 


As discussed in the introduction of this article, modern image-based 
characterization techniques enable measuring grains-size distributions 
of polycrystalline materials empirically. Although analytically given 
grains-size distributions have their merits, relying upon empirically 
given GSDs permits a close comparison to real-world specimens, in 
particular concerning their inherent imperfections and deviations from 
theoretical GSDs. 

To demonstrate the capabilities of the proposed methodology, we use 
the emprical GSDs determined by Dobrich et al. (2004) as our target 
GSDs for generating synthetic microstructures. More precisely, we use 
histogram bins for describing the desired GSD, see Fig. 3.8. 
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Figure 3.8: Comparison of targeted and reached bin frequencies for the GSD of Dobrich 
et al. (2004). (a) 250 grains, (b) 500 grains, (c) 1000 grains. 
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Different centroidal periodic Laguerre tessellations were generated with 
an increasing number of grains using Barzilai-Borwein in combination 
with the Anderson-Lloyd method. The resulting microstructures are 
shown in Fig. 3.9, and the run-times are collected in Tab. 3.12. We see 
that even the largest structure with 8000 grains required only slightly 
more than two minutes to compute. 


(b) (c) 


(d) (e) (f) 
Figure 3.9: Generated microstructures with empirical GSD of Dobrich et al. (2004). 


(a) 250 grains, (b) 500 grains, (c) 1000 grains, (d) 2000 grains, (e) 4000 grains, 
(f) 8000 grains. 
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Number of grains | 250 500 1000 2000 4000 8000 
run-time in s 0.53 1.62 3.04 10.57 33.15 133.51 


Table 3.12: Run-times for different numbers of grains and the empirical GSDs of Débrich 
et al. (2004) using BB and Alg. 2. 


Of particular importance to us is the ability of the methodology to match 
any prescribed empirical GSD. For this reason, we take a look at Fig. 3.8, 
where both the target and the realized histograms are shown for 250, 
500 and 1000 grains. For 250 and 500 grains, differences between the 
histograms are visible, although on a level probably within engineering 
accuracy. For 1000 grains, the target and the realized histograms are 
difficult to distinguish by eye. 

To gain more insight, we plotted the absolute errors (in %) between 
the desired and targeted frequencies per bin, cf Fig. 3.10. We see a 
maximum error of about 0.5% for 250 grains, 0.3% for 500 grains and 
below 0.15% for 1000 grains. For a higher grain count, the error decreases 
further. More precisely, we see that the error decreases as 1/K, where 
K denotes the number of grains, i.e., the error reduces roughly by 50% 
when doubling the number of grains. This is no surprise, but a result 
of our quasirandom sampling of the grain sizes, which theoretically 
predicts a 1/K decrease of error upon sampling (Asmussen and Glynn, 
2007). In contrast, a purely random sampling only reduces the expected 
error by a factor of 1/ VK. To gain some intuition into this difference, 
essentially to arrive at the same error for random sampling that we 
obtained with quasirandom sampling and 1000 grains would require 
10007, i.e., one million, grains, on average. Last but not least let us 
remark that the maximum error decreases linearly, but the bin-wise error 
might not. Indeed, taking a look at the second-smallest bin, the error 
only decreases for 2000 grains compared to 250 grains. However, the 
corresponding error is already small for 250 grains. 


78 


3.5 Examples with a higher level of complexity 


0.5 4 0.5 F 4 
0.4 =| 0.4 F | 
0.3 + 0.3 F 3 


absolute frequency error in % 
absolute frequency error in % 


0.8 14 16 0.2 04 06 0.8 1 12 14 16 
equivalent normalized radius equivalent normalized radius 
(a) (b) 
T T T T 
x 0.5 F J xæ 0.5 H + 
8 04H 7 8 047 5 
a a 
© © 
> > 
S 03h 48034 | 
S S 
3 A 
a 0.24 4 È 02} 4 
£ g 
2 01H 4% 01H J 
a a 
0 = = 07 n 
0.2 04 06 0.8 1 12 14 16 02 04 06 0.8 1 12 14 16 
equivalent normalized radius equivalent normalized radius 
(c) (d) 
T T T T T T T T T T T 
x 0.5 F Jæ 051 = 
5 Aal {15 04, | 
8 = 
© © 
5 > 
= 03f + 3 0.3 -| 
g Sg 
a a 
D D 
8 02F 4 & 0.24 | 
2 2 
2 016 | 2 01f 4 
a a 
0 m - Q m a | 
0.2 04 06 0.8 1 12 14 16 0.2 04 06 0.8 1 12 14 16 
equivalent normalized radius equivalent normalized radius 
(e) (£) 


Figure 3.10: Absolute frequency error of the microstructure realizations of Fig. 3.9, as 
a result of the quasi-random sampling, compared to the target values of Döbrich et al. 
(2004). (a) 250 grains, (b) 500 grains, (c) 1000 grains, (d) 2000 grains, (e) 4000 grains, 
(f) 8000 grains. 
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3.6 Conclusion 


This work has been devoted to generating geometry-based microstruc- 
ture models for polycrystalline materials in terms of Laguerre tessella- 
tions. According to Ball and Carstensen (2015) Laguerre tessellations 
describe polycrystalline microstructures with convex grains. Put differ- 
ently, in order to describe a geometrically more complex polycrystalline 
microstructure, non-convex grains have to be considered. 

In a recent work, Bourne et al. (2020) have shown that, for any given set 
of seed points and prescribed volume fractions, a corresponding set of 
Laguerre weights realizing a Laguerre tessellation with the prescribed 
volume fractions may be determined by solving a convex program. As 
a consequence, using block-box optimization methods, Bourne et al. 
(2020) were able to generate polycrystalline microstructures of industrial 
complexity in the order of minutes. 

In this article, we retained the basic algorithmic framework of Bourne 
et al. (2020) for computing the weights of a periodic Laguerre tessellation 
for prescribed seed locations and volume fractions, and exploited the 
benefits of modern non-monotone methods of gradient type, i.e., the 
Malitsky-Mishchenko method and the Barzilai-Borwein scheme, in the 
context of semi-discrete optimal transport. 

In the numerical experiments we saw that, for small to moderate grain 
count and low standard deviation of the grain-size distribution, New- 
ton’s method outperforms all other methods considered. For larger 
number of grains or for larger spread of the grain sizes, the Barzilai- 
Borwein method proved superior in terms of run-time. As the run-time 
for small- to medium-sized problems is small anyway, the Barzilai- 
Borwein method may serve as a general-purpose method for computing 
the weights of Laguerre tessellations. 

To increase the physical realism, we also contributed to computing cen- 
troidal Laguerre diagrams of prescribed volume fractions by suggesting 
to combine the classical Lloyd’s algorithm to Anderson acceleration, a 
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general-purpose acceleration scheme for fixed-point methods. For all 
problems considered, the Anderson-accelerated Lloyd’s method proved 
superior to the classical Lloyd’s method. Anderson acceleration might 
also prove more powerful if higher accuracy for centering is required, 
in particular when compared to BFGS-type centering (Liu et al., 2016; 
Hateley et al., 2015). However, in the context of polycrystalline mi- 
crostructures, such a high level of accuracy appears not to be necessary. 
Allin all, combining the improved convex optimization solver, the cen- 
tering and an OpenMP-parallel implementation permits generating poly- 
crystalline microstructures of industrial complexity and relevance in a 
matter of seconds to minutes, at most, reducing the run-time and increas- 
ing the number of grains compared to the work of Bourne et al. (2020). 

In particular, for grain-size distribution with large standard deviation 
of the grain sizes, it appears imperative to accept iterates with van- 
ishing Laguerre cells. Classically, backtracking methods are used to 
avoid empty Laguerre cells during the iterative process. The proposed 
gradient methods avoid such backtracking - in fact, their performance 
decreases when backtracking is used. Furthermore, we also proposed 
a simple regularization strategy for Newton’s method in case of van- 
ishing Laguerre cell volume which leads to a robust and competitive 
computational scheme. 

Our solution schemes may be further improved by using multi-level 
strategies, as pioneered by Mérigot (2011) and Lévy (2015). However, 
for applications to micromechanics, the number of grains to consider 
typically does not exceed a few thousand. 

Notice that we describe the morphological texture of a polycrystalline 
aggregate in terms of the seed locations and the volume fractions, taking 
into account only one-point statistics (Torquato, 2002). We have shown, 
by example, that microstructures with a spatial gradient in the grain-size 
distribution may be generated. Also, experimentally determined grain- 
size distributions may be reproduced with high accuracy, provided the 
equivalent radii are drawn in a pseudo-random fashion. Still, it might 
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be interesting to see whether the algorithms could be extended to more 
general morphological features, e.g., sphericity. 

Last but not least it is clear that, in order to conduct computational 
experiments, the generated polycrystalline microstructures need to be 
furnished with crystal orientations per grain, leading to possible crys- 
tallographic texture of the material under consideration. However, the 
latter is typically realized by a post-processing procedure, see, e.g., Quey 
et al. (2018), and is thus complementary to the approach presented. 


82 


Chapter 4 


Generating polycrystalline 
microstructures with prescribed 
texture coefficients ' 


4.1 Introduction 


The methods proposed in Chap. 3 enable the generation of polycrys- 
talline microstructures of industrial complexity within several seconds 
up to a few minutes. As outlined in Sec. 2.4 the next step in the gener- 
ation of synthetic polycrystalline microstructure is assigning crystallo- 
graphic orientations to each grain. 

In this chapter, we propose a method based on condensing the informa- 
tion carried by the crystallite orientation distribution function (CODF) 
via the coefficients of a tensorial series expansion (Guidi et al., 1970; 
Adams et al., 1992). As these tensorial coefficients are easy to compute 
from experimental data, directly linked to the bounds of mechanical 
properties (Böhlke and Bertram, 2003; Bohlke and Lobos, 2014) and easily 
applied in microstructural sensitive design (Fullwood et al., 2010; Lobos 
and Bohlke, 2015), we prefer them to an approach involving spherical 
harmonics (Roe, 1965; Bunge, 1982). Moreover, Bohl ke (2005) introduced 


! This chapter is based on the publication by Kuhn et al. (2022) whereas minor 
typographical and formatting changes have been made for cohesion of the manuscript 
and the readers convenience 
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a method for estimating the crystallite orientation distribution function 
for a finite number of given texture coefficients and Junk et al. (2012) 
analyzed this approach in the context of maximum entropy moment 
problems. Böhlke (2006) derived a hierarchy of evolution equations for 
the texture coefficients under the Taylor-Voigt assumption of vanishing 
strain fluctuations on the microscale. Motivated by these observations, 
we propose to use a limited set of tensorial texture coefficients as the 
input of a microstructure generator, with the goal to obtain crystallo- 
graphic orientations for each crystallite matching the prescribed texture 
coefficients. We follow a two-step procedure. First, the orientations are 
sampled randomly. In a second step, these orientations are corrected 
in terms of a gradient-based optimization technique, which ensures 
that the resulting texture coefficients match the prescribed ones up to 
a given tolerance. To introduce the Texture coefficient Optimization for 
Prescribing orientations (TOP) method, we briefly revisit the notation of 
rotations in Sec. 4.2.1. In Sec. 4.2.2, we outline the problem formulation 
and the solution scheme. We investigate the capabilities of the proposed 
method, by comparing it to state-of-the-art algorithms in Sec. 4.3. 


4.2 Background on modeling 
and optimization 


4.2.1 Representing the texture 


Describing the deformation behavior of a polycrystal, i.e., an agglomer- 
ate of crystals, by crystal plasticity requires taking the distinct orientation 
of each crystallite into account. For instance the stiffness tensor C in 
Hooke’s law (2.31) depends on the crystal orientation. To describe the 
orientations in a polycrystalline material, for each crystallite we use a 
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proper orthogonal tensor 
3 
Q= 9% 8 &, (4.1) 
i=1 


where (€1, €2, 3) and (91, 92, 93) represent the fixed orthonormal basis 
of the sample and the crystallite, respectively. Thus, the orientation 
of a crystallite is encoded by the rotation from the crystal coordinate 
system into the sample coordinate system. All orientation tensors Q 
are elements of the group of proper rotations in three dimensions, i.e., 
Q e SO(3). In the following, we use the expressions rotation and 
orientation interchangeably. 

For a given polycrystal, the orientation may be succinctly described in 
terms of the crystallite orientation distribution function (CODF) f, a 
probability distribution on SO(3) w.r.t. the Haar Measure dQ (normal- 
ized to unity). More precisely, the CODF f is non-negative, 


f(Q) 20 Y QeSO(3), (4.2) 


and normalized 


/ (9) dQ =1. (43) 
so(3) 


For a (measurable) subset S C SO(3) of orientations, the expression 
Js F(Q) dQ computes the probability to find orientations contained in 
the set S. Due to the invariance properties of the Haar measure (Gel’fand 
et al., 2018), the invariance property 


/ f(Q)4Q-= | f(QQo)dQ holds for all Qo € SO(3). 
S0(3) SO(3) üa 
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Moreover, the CODF reflects the underlying symmetries of the crystals 
forming the aggregate, i.e., 


f(Q) = f(QH®) holds forall H? € S°CSO(3), (4.5) 


where S7 denotes the (discrete) symmetry group of the crystals, as the 
orientation states Q and QH“ correspond to the same physical orienta- 
tion state of the crystallite. As a result of from a forming process, the 
sample itself may possess a certain symmetry, encoded by a symmetry 
group S5. This is reflected by the CODF in terms of the condition 


f(Q) = f(H°Q) V H? € S° c 80/3), (4.6) 


where $° denotes the symmetry group of the sample. For the sake 
of readability, we will only consider the case of cubic crystals and 
triclinic sample symmetry in the following (Morawiec, 2003, Ch. 3). Our 
approach permits a straightforward extension to the general case with 
arbitrary crystal and sample symmetries, see Zheng and Zou (2001a;b). 
If the CODE is not uniform, i.e., f(Q) # 1 for some Q € SO(3), then the 
material will be said to possess a crystallographic texture. 

Working with the full CODF is oftentimes impractical. Guidi et al. (1970) 
proposed a way to condense the encoded information. More precisely, 
any square integrable function f may be expressed in terms of a tensorial 
Fourier series (Guidi et al., 1970; Adams et al., 1992) 


f(Q)=14 90 Via, Fia) (Q); (4.7) 
i=1 


involving the tensorial Fourier, or texture, coefficients Via) and the 
tensor functions Fo) (Böhlke, 2005; Fernández and Böhlke, 2019). The 
subscript i denotes the number of linear independent harmonic cubic 
tensors of rank a. Harmonic tensors, denoted by a prime (-)’, are 
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completely symmetric and traceless, i.e., the relations 


3 
Ai ikl Aiki Ajik “Ty So Abit =0 (4.8) 
i=1 


reduce the number of degrees of freedom to 2a + 1 (Böhlke, 2005). 


The tensor functions Bias (Q) in equation (4.7) arise by rotating suitable 


reference tensors Das (Guidi et al., 1970) 


F(a.) (Q) = Q * Tha) (4.9) 
where 
QxT = Tij.. p(Qe:i) 8 (Qe;) 8 ... (Qep) (4.10) 


denotes the Rayleigh product, i.e., the rotation of a tensor by Q. Without 
loss of generality and following Böhlke et al. (2010) and Dyck and Böhlke 
(2020), we normalize the Frobenius norm of the reference tensors to 


Eal =2a+1 (4.11) 


instead of IT, | = 1 (as done by Guidi et al. (1970)). Interpreting 
the texture coefficients Via) as a convex combination of normalized 


AR 


reference tensors of single crystal states (Fernändez and Böhlke, 2019), 
we compute the tensorial texture coefficients in the case of discrete 
(experimental) orientation data as (Böhlke, 2006) 


K 
1 
Vian (Qr Qr) = 5, 7 NO ce Qe * Tas (4.12) 
: k=1 


where cp denotes the volume fraction of orientation Q, among K orien- 
tations. From equations (4.11) and (4.12) it follows that the Frobenius 
norm of all texture coefficients V}, lies within the interval [0, 1]. Böhlke 
(2006) and Böhlke and Lobos (2014) use the norm to measure the degree 
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of anisotropy. For a completely uniform CODE, all texture coefficients 
vanish, whereas for the case of single crystals the norm of all texture 
coefficients is equal to one. 

As we seek a compact representation of the CODF in the following, 
we will restrict to the texture coefficients up to rank six. As the cubic 
reference tensor of rank two is zero and, because of the cubic crystal 
symmetry, odd-rank reference tensors up to a rank of eight vanish (Guidi 
et al., 1970), we focus on the texture coefficients of rank four and six. For 
a precise overview and thorough discussion of texture coefficients in a 
more general context, the reader is referred to the work of Fernandez 
and Böhlke (2019). 


4.2.2 Texture coefficient optimization 
for prescribing orientations (TOP) 


To create digital representations of polycrystalline microstructures, it 
is not sufficient to solely match the grain morphology. In addition, we 
have to take the orientation state, i.e., the CODF, into account (Kocks 
et al., 2000). Many tools either rely on simple model CODFs to generate 
orientations (Quey et al., 2011; 2018), need (possibly vast) experimental 
data (Groeber and Jackson, 2014) or at least a representation of the 
complete CODF (Eisenlohr and Roters, 2008) to generate orientations for 
digital microstructures. In this section, we propose a method to generate 
orientations based on tensorial texture coefficients. For a given unit cell, 
subdivided into individual grains, our goal is to prescribe the orientation 
per grain in such a way that the resulting texture coefficients V/) of 
the unit cell match the prescribed ones V/,,) up to a given tolerance tol, 
thus approximating the underlying CODF. To this end, we formulate our 
objective function as the difference in independent components between 
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the current and the desired texture coefficients 


Omax 


(Qe) = 2 IV, — Vig Qe) IP, (4.13) 


where Qc = (Qı,..:,Qx) is a vector of orientations (one for each of 
the K cells) and amax denotes the maximum rank of the considered 
texture coefficients. 

Starting from a randomly initialized set of orientations, we seek to 
minimize the objective function £(Qc). The objective function Zis contin- 
uously differentiable, and it natural to use gradient-based optimization 
techniques. Indeed, the objective function ¢ is actually a polynomial in 
the components of the individual grain orientations @,,. In particular, 
the function £ is infinitely often differentiable. 

The simplest conceivable approach proceeds via gradient descent, i.e., 
following the direction of steepest descent. Please keep in mind that 
the objective function £ is defined (4.13) on the K-fold product of SO(3). 
The space of special orthogonal tensors SO(3) forms a non-linear subset 
of the space of second-order tensors. In particular, a simple explicit 
Euler discretization of gradient descent does not work, as the next trial 
point does not necessarily lead to vector whose components satisfy the 
constraints of SO(3) (Taylor and Kriegman, 1994; Makinen, 2008). 

This is illustrated in Fig. 4.la. On a curved space with Riemannian 
metric, the natural extension of straight lines are so-called geodesics, 
which emanate from a point in a specific (tangent) direction by par- 
allel translation. On a (compact, matrix) Lie group with its natural 
Riemannian metric (the Killing form), following the geodesics may 
be computed in terms of the matrix exponential. In case of SO(3), 
this reduces to Rodrigues’ formula 


(4.14) 
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describing a rotation around an axis u by an angle 0 and where we set 
w = ĝu as well as 


0 wg wa 
J(w) = w3 0 -w |. (4.15) 
—W2 Wy 0 


The gradient descent scheme works as follows, where a step size t > 0 is 
fixed beforehand. Suppose the i-th iterate Qi, = (Q4, .. . , Qi) is given 
(i =0,1,...). Then, we investigate the function 


£;(w) = &(Qiexp (J(wi)),..., Qxexp (J(wx))), (4.16) 


where w = (w1, ... wg) € R°*. Computing VL;(0) € R° by the chain 
rule, gradient descent proceeds via 


a = (Qiexp (J(-t[V4i(0)]1)),---, Qexp (J(—t[V4i(0)]))) 
(4.17) 


(a) (b) 


Figure 4.1: (a) Gradient descent on a linear space vs. (b) descent along a geodesic (dashed 
line) on the manifold SO(3). 
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4.3 Computational investigations 


4.3.1 Setup 


In the following sections, we wish to provide insights into the per- 
formance of the proposed TOP method. As a first step, we consider 
linear elastic material behavior, i.e., we investigate the effective stiffness. 
Secondly, we study cyclic stress-strain hysteresis. 

To create the morphology of the microstructures under investigation, 
we use the algorithm described in Kuhn et al. (2020) to generate digital 
polycrystalline microstructures with prescribed volume fractions. For 
the morphology we consider two cases, a unique and log-normal grain 
size distribution (GSD). The former means that all grains have the same 
volume, i.e., V} = 1/a, where G denotes the total number of grains 
in the volume element. Restricting to a unique grain size permits us 
to study the influence of the individual grain orientations exclusively. 
Due to their frequent occurrence in experiments (Spettl et al., 2014), we 
also investigate microstructures with an equivalent diameter following 
a log-normal grain size distribution with mean equal to unity and a 
standard deviation of 0.15, see Kuhn et al. (2020). 

For a fixed morphology, we furnish the grains with orientations, where 
we investigate three different CODFs, namely a uniform, one with a 
slight texture and one with an increased texture. For the former, we 
compare the accuracy of TOP and the algorithm proposed by Quey et al. 
(2011; 2018), integrated into the polycrystal generation software Neper. 
In addition, we include a random orientation sampling (realized using 
the scipy implementation for sampling the Haar distribution (Stewart, 
1980)) as a benchmark. For the textured CODFs, we compare the TOP 
method to random sampling from discrete orientation measurements, 
which is a common practice (Gillner and Münstermann, 2017; Tu et al., 
2019). For the textured CODFs and a log-normal GSD, we additionally 
consider the Texture Discretization Technique (TDT) algorithm proposed 
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by Melchior and Delannay (2006), which, in a first step, samples ori- 
entations using the method proposed by Tóth and Van Houtte (1970). 
In a subsequent clustering step, a binary look-up table is computed by 
evaluating the misorientation of each pair of sampled orientations. If 
this misorientation is below a chosen threshold value, the corresponding 
entry in the look-up table is set to 1 otherwise to 0. To assign orientations 
to grains, each grain is associated with a number of so-called elementary 
volumes according to their size, which is used to find orientations in 
the look-up table with at least this number of orientations having low 
misorientation. The corresponding crystallographic grain orientation is 
then the average of orientations with low misorientation to each other. 
The parameters, i.e., the number of elementary volumes per grain and 
the threshold value for the misorientation, have to be chosen judiciously. 
The TOP method is implemented in Python with Cython extension 
following the optimization procedure outlined in section 4.2.2. Unless 
otherwise specified, we use a tolerance of tol = 1078 to solve the 
optimization problem and consider texture coefficients up to rank six. 
The material model described in Sec. 2.2.2 is implemented in a user- 
material-subroutine (UMAT). The coefficients of the elastic stiffness 
tensor are taken from the literature (Wu et al., 2017; Schäfer et al., 2019a), 
whereas the critical resolved shear stress, assumed to be identical for 
all slip systems, and the parameters of the kinematic hardening model 
were fitted to experimental stress-strain hysteresis of the steel C45 using 
Bayesian optimization (Kuhn et al., 2021). The complete set of used 
model parameters is summarized in Tab. 4.1. 
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cubic stiffness (in Voigt notation) Cy, = 231.4GPa C2 = 134.7 GPa C44 = 116.4 GPa 


flow rule 40 = 0.0015”! m = 100 
Ohno-Wang model parameters M=8 Te = 81.4 GPa 
A= 6910.0 MPa B = 6734.0 
lattice type BCC 
slip systems {110}(111) 
Table 4.1: Parameters used for the crystal plasticity 


model (Wu et al., 2017; Schäfer et al., 2019a). 


To efficiently compute the effective stiffness as well as the macroscopic 
stress-strain hysteresis, we use the FFT-based solver FeelMath (Fraun- 
hofer ITWM, 2021; Wicht et al., 2020a;b). For the stiffness computa- 
tions we rely on the conjugate gradient method (Zeman et al., 2010; 
Brisard and Dormieux, 2010), whereas for the non-linear problem we 
use a Newton-CG method (Gélébart and Mondon-Cancel, 2013; Ka- 
bel et al., 2014). For both problems we use the Moulinec-Suquet dis- 
cretization (Moulinec and Suquet, 1994; 1998). For a perspective of 
solution schemes and discretizations, we refer to the recent review article 
by Schneider (2021). By default, we carry out the computations on peri- 
odic microstructures, discretized by 64° voxels. Please note that we apply 
periodic boundary conditions to compute the stiffness and hysteresis. 


4.3.2 Linear elastic stiffness 


In this section, we study the effect of different orientation-sampling 
techniques on the effective stiffness of polycrystalline microstructures. 
In order to minimize the influence of the underlying microstructure 
morphology, we use a fixed grain microstructure for each realization and 
all orientation sampling methods. This is illustrated in Fig. 4.2, where 
we show the results of different sampling techniques for a fixed grain 
structure with grains of identical volume. 
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u 


10017 on] 


(a) (b) 


Figure 4.2: Furnishing a grain microstructure with 128 grains of identical volume with 
orientations sampled by three different methods. The colors correspond to the inverse- 
pole-figure color key in 010 direction (Nolze and Hielscher, 2016). (a) TOP, (b) Neper, 
(c) Random. 


Uniform CODF We start with the case of a unique grain-size and a uni- 
form orientation distribution, corresponding to mechanically isotropic 
behavior (Krawietz, 1999; Bertram et al., 2000; Böhlke and Bertram, 2001). 
For the results to be representative, it is necessary to determine the 
number of grains which ensure an isotropic effective material response, 
see for example Kanit et al. (2003) and Yang et al. (2019). In this spirit, 
we investigate microstructures with an increasing number of grains and 
study their effective stiffness. 

As discussed in Sec. 4.2.1, for a uniform CODE, all texture coefficients 
vanish, i.e., 

Vig) =0 (4.18) 
holds for all considered texture coefficients. To quantify the anisotropy 
of the stiffness tensor we compare to the best approximation by an 
isotropic tensor (see equation (4.22)), i.e., we project the computed 
stiffness tensor onto the space of isotropic tensors of fourth order. For a 
detailed discussion see the work by Fedorov (1968) and Arts (1993). 
We compute the mean stiffness 


1 N 
Ce = N 5 Cen (4.19) 
n=1 
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of N = 10 realizations and extract the Lamé constants via 


pP = (Co,a4 a Co,55 ar Ca,66) (4.20) 


Ol = w| rm 


NPP = = (Ca 12 + Ca,13 + Ca,23 + Ca,21 + Ca,31 + Ca,32), (4.21) 
where Cg,i; denotes the ij-th entry of the stiffness tensor in Voigt 
notation (Cavallini, 1999). 
Using the best isotropic approximation C's (PP, APP) based on the 
extracted Lamé constants, we introduce the isotropy error &® via 

_ [Ce (utPP, APP) — Call 


S (Ca, Ca, N) = ICa] ; (4.22) 
e 


measuring the degree of anisotropy present in the computed stiffness. 
For an increasing number of grains 


G € 132, 64, 96, 128, 256, 384, 512, 768, 1024, 1536}, 


we show the resulting isotropy error for the three different orientation 
sampling methods in Fig. 4.3a. We observe a decreasing isotropy error 
for all methods with an increasing number of considered grains. All 
methods decrease the isotropy error at a similar rate. However, they 
differ in the initial error level. For instance, all sampling techniques reach 
a low isotropy error for 1536 grains, namely 0.251%, 0.041% and 0.027% 
for random sampling, the Neper and TOP methods, respectively. To 
reach a mean error below 1%, the microstructure has to consist of more 
than 64 grains if the orientations are sampled randomly or generated 
by Neper. For all investigated grain counts, the TOP method produces 
the lowest isotropy error. Neper starts with a substantially higher error 
(by roughly one order of magnitude) at low grain counts and reaches 
a similar performance to TOP for more than 300 grains. For the naive 
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random sampling, the isotropy error has a quite large offset to the more 
involved algorithms. 


— Neper — TOP — Random 
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Figure 4.3: Isotropy and total error for effective stiffnesses computed from microstructures 
with uniform orientations and unique grain size distribution. (a) Isotropy error, 
(b) Total error. 


In addition to evaluating the degree of isotropy, we investigate the 
deviation from the effective, infinite-volume stiffness. As our ground 
truth, we consider the mean of ten apparent stiffnesses, each computed 
using volume elements consisting of 10000 grains and discretized by 
128° voxels (see Fig. 4.4a). 
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u 
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(b) 


Figure 4.4: One realization of microstructures with with TOP based orientations. The color 
corresponds to the ipf color key in 100 direction. (a) 10 000 grains, (b) 1024 grains. 


The orientations are sampled using the TOP method. The mean 
stiffness and the 95% confidence intervals, computed via Student’s 
t-distribution (Student, 1908), are given in Tab. 4.2. The isotropy error of 
the mean stiffness is 6° = 0.012%. 


276.90 +0.07 112.00 + 0.08 112.00 + 0.09 0.00 0.10 0.00 0.08 0.00 0.08 
— 276.90 + 0.07 112.00 0.06 0.00 0.08 0.00 0.05 0.00 + 0.08 
= _ 276.90 + 0.10 0.00 + 0.07 0.00 + 0.05 0.00 + 0.05 
= = _ 82.50 + 0.15 0.00 + 0.06 0.00 + 0.08 
— = = = 82.50 + 0.06 0.00 + 0.07 
= = = = 82.50 + 0.11 


Table 4.2: Mean and 95% confidence intervals in GPa for the stiffness in Voigt’snotation 
computed by averaging ten realizations of microstructures with 10000 grains and 
uniformly distributed TOP orientations. 


We define the total error tt as the mean relative error between the 
stiffness of each realization and the one given in Tab. 4.2, i.e., 


ste IC Conlı (4.23) 
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For the total error, shown in Fig. 4.3b, we make similar observations as 
for the isotropy error. All of the methods decrease the total error at a 
similar rate, but differ initially. For the case of orientations generated by 
Neper, the error for 32 grains is ö'° = 1.17% and therefore roughly twice 
as large compared to the TOP method with 6'°! = 0.521%. To reach a 
similar error with randomly sampled orientations about 768 grains have 
to be considered. Random sampling leads to a mean error of 0.677% for 
1536 grains. The errors produced by the Neper and TOP method are 
similar to each other with 0.073% and 0.078%, respectively. 

To understand the similar rates of error decrease more thoroughly, it is 
helpful to decompose the total error ôt into two contributions (Kanit 
et al., 2003). The first part is the random error and quantifies the 
inaccuracy associated with working on a reduced representation of the 
ground truth. The second contribution quantifies artificial long-range 
correlations introduced by working on periodic microstructures (Kanit 
et al., 2003; Schneider et al., 2021). We attribute the visible offset in 
Figs. 4.3 and 4.6 to the random error, as we use the same geometric 
representations for each orientation sampling method. Thus, a smaller 
random error is achieved by the TOP method and further reduction of 
the total error 5'* is attributed to increasing the cell-size, i.e., increasing 
the number of grains. 

Up to this point, our investigations were based on a polycrystal with 
a unique GSD, i.e., a unique grain size. However, to account for the 
influence of the grain size on the mechanical response, it may often be 
necessary, and therefore desirable, to match more realistic grain size 
distributions when generating synthetic polycrystalline microstructures. 
Thus, we turn to polycrystals with a log-normal GSD, as typically ob- 
served in real-world samples (Spettl et al., 2014), with a mean equivalent 
diameter equal to unity and a standard deviation of 0.15. Fig. 4.5 shows 
an example of a microstructure consisting of 128 grains, equipped with 
orientations from the three different sampling techniques. 
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Figure 4.5: Polycrystalline microstructure realizations with 128 grains following a log- 
normal grain size distribution and orientations generated by three different methods. The 
colors correspond to the ipf color key in 010 direction. (a) TOP, (b) Neper, (c) Random. 


The isotropy as well as the total error for different sampling methods are 
shown in Fig. 4.6. For the log-normal GSD, we register a notable decrease 
in accuracy for the Neper sampling method compared to the unique 
grain size. Considering the case of 32 grains, the isotropy error and total 
error increase from 1.22% and 1.17% to 1.62% and 2.19%, respectively. 
This loss of accuracy persists, even for larger grain counts, e.g., for 
1536 the total error for the unique and log-normal GSD are 0.07% and 
0.37%, respectively. For randomly sampling the Haar distribution, the 
influence of a log-normal grain size distribution is smaller. For instance 
the biggest difference in the isotropy error is 1.23% and 1.60% for the 
unique and log-normal GSD and 32 grains, respectively. The proposed 
TOP method takes the volume fraction of each grain into account in an 
explicit way when optimizing the orientations. This results in strikingly 
similar error levels for both the unique and the log-normal case. Whereas 
the total error values realized by microstructures with 32 grains differ 
slightly for the unique and log-normal case, the resulting isotropy error 
is 6° = 0.16% for both GSDs. 
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Figure 4.6: Isotropy and total error for the effective stiffnesses computed from microstruc- 
tures with uniform orientations and log-normal grain size distribution. (a) Isotropy error, 
(b) Total error. 


We investigate the influence of the maximum rank of the texture coeffi- 
cients considered in our optimization scheme in the case of a uniform 
orientation distribution. For this purpose, we consider the case of ten 
microstructures consisting of 1024 grains, see Fig. 4.4b for an example of 
one realization. We study two different cases: Using solely the texture 
coefficient of rank four to optimize the orientations and considering the 
coefficients of rank four and six. For these two cases, the isotropy and 
total error are shown in Fig. 4.7 for different values of the tolerance tol 
used in the optimization procedure. Both errors show a decrease up to 
a value of tol = 1076, after which the resulting errors do not change. 
Interestingly, considering only the texture coefficient of rank 4 appears 
to be beneficial. As all total errors are below 0.1% and all isotropy errors 
even below 0.03%, using texture coefficients of rank four and solving 
the problem up to a tolerance of tol = 10~° is sufficient for the case of 
a uniform orientation distribution when considering 1024 grains and 
solely the macroscopic stiffness is of interest. 
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Figure 4.7: Comparison of the 5'° and ô% for the stiffness computed for 1024 grains with 
TOP orientations and different considered texture coefficients plotted on a finely resolved 
y-axis. (a) Isotropy error, (b) Total error. 


Textured CODF To further investigate the capabilities of the TOP 
method, we turn to a non-uniform CODE i.e., a textured polycrystal. 
The prescribed CODF was generated by MTex (Hielscher and Schaeben, 
2008), taken from the MTex documentation (MTex, 2017), see Fig. 4.8 
for the corresponding pole figures. As MTex allows the sampling of 
CODEFs, we draw 50 000 samples at random for computing the texture 
coefficients, assuming the same weight for each sample. 

As a ground truth we define the mean stiffness of ten realizations, 
each with 10000 grains. The resulting stiffness for TOP orientations 
is given, with its respective 95% confidence intervals, in Tab. 4.3. The 
isotropy error of this stiffness computes to 5° = 5.54%, i.e., a slight 
anisotropy appears. 
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Figure 4.8: Pole figures of the (generated) textured CODF (MTex, 2017). 


269.70 + 0.07 117.10 +0.06 113.90 + 0.05 0.00 + 0.07 —2.60 + 0.11 0.10 + 0.07 
= 268.60 + 0.08 115.00 + 0.06 0.10+0.08 —0.10 +0.06 —0.10 + 0.05 
Z— _ 271.90 + 0.07 —0.10 0.09 2.70 0.06 0.00 0.07 
= = = 85.80 + 0.14 0.00 +0.10 —0.10 + 0.08 
= = 84.60 + 0.16 0.00 0.05 
= = = = = 88.10 + 0.16 


Table 4.3: Mean and 95% confidence intervals in GPa for the stiffness in Voigt notation 
computed by averaging ten realizations of microstructures with 10000 grains and TOP 
orientations for a synthetic CODF. 


For this texture, we investigate the approximation quality of the stiffness 
for a varying number of grains, each with identical volume. The total 
error for randomly sampling from the generated orientations and using 
texture coefficients is shown in Fig. 4.94. Randomly sampling the 
generated orientations gives a total error of ĝt = 4.47% for 32 grains 
and reaches an error below 1% for 936 orientations. The mean total error 
achieved by the TOP method of 6'°t = 0.51% for 32 grains actually lies 
below the error value achieved by randomly sampling 1536 orientations. 
For the latter number of grains, TOP achieves an error 5° = 0.07%. 
This difference is attributed to the notable offset between the random 
sampling and TOP method, as both decrease ôt* with the same rate. 

Let us consider the case of a log-normal grain size distribution. For the 
TDT algorithm we assign eight elementary volumes to the smallest grain 
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and increase the number of elementary volumes for each grain according 
to its size. We set the threshold misorientation value to 5°. In Fig. 4.9, we 
provide the total error 5‘ for ten realizations with a varying number of 
grains. For the case of randomly sampling from given orientations, we 
observe a slight increase in the error value induced by the underlying 
log-normal grain size distribution. For instance, for a microstructures 
with 32 grains, the mean error is 1% = 3.12% and 6" = 3.74% for the 
unique and log-normal GSD, respectively. This effect decreases when a 
larger number of grains is considered, as the effect of a single, large grain 
with specific orientation on the overall response decreases. In contrast, 
the TOP method is not adversely affected. Indeed, during optimization, 
the volume fraction is explicitly taken into account when computing 
the texture coefficients, see equation (4.12). Using the TDT algorithm 
results in a lower total error than random sampling for all grain numbers 
considered. For instance, for a 64-grain microstructure, the total error 
is 61 = 3.75% and ö'%* = 2.77% for random sampling and the TDT, 
respectively. For both algorithms the error decreases with a similar rate 
as for the proposed TOP algorithm, whereas they both result in higher 
total errors than using TOP. 
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Figure 4.9: Total error 5'*' of the effective stiffness computed for microstructures with a 
unique and log-normal grain size distribution and a synthetic CODF. (a) Unique grain size 
distribution, (b) Log-normal grain size distribution. 


Highly textured CODF In practical applications, e.g., cold rolled steel, 
the intensities in the pole figure may reach values as high as ten. To 
investigate this scenario, we next consider a case with an increased 
texture in the CODF. We rely on synthetically generating a CODF us- 
ing MTex (Hielscher and Schaeben, 2008) and show the resulting pole 
figures in Fig. 4.10. 
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Figure 4.10: Pole figures of the (generated) CODF with increased texture (MTex, 2017). 


For the ground truth we proceed in the same way as for the slightly 
textured CODF, using ten microstructures with 10 000 grains of equal vol- 
ume, equipped with orientations from the TOP method to compute the 
mean stiffness. For this case, the mean and the 95% confidence intervals 
are given in Tab. 4.4. The isotropy error is ö'°° = 18.78% which is more 
than three times the error of the slightly textured case, i.e., 5° = 5.54%. 


257.10 +0.06 115.40 +0.07 128.30 + 0.04 0.00 + 0.06 —0.00 + 0.03 —0.20 + 0.03 
_ 260.60 + 0.07 124.80 + 0.04 0.00 + 0.04 0.00 + 0.03 0.20 + 0.05 
_ 247.60 + 0.05  —0.00 + 0.05 0.00+ 0.04  —0.00 + 0.04 

= 98.40 + 0.09 —0.10 + 0.06 0.00 + 0.10 

— — — _ 103.70 + 0.09  —0.00 + 0.07 
= 85.60 + 0.17 


Table 4.4: Mean and 95% confidence intervals in GPa for the stiffness in Voigt notation 
computed by averaging ten realizations of microstructures with 10000 grains and TOP 
orientations for a synthetic CODF. 


First, we investigate the case of a unique grain size distribution and 
show the resulting total error in Fig. 4.1la. For TOP and random 
sampling, the error decreases with a similar rate, which is consistent 
with our observations in the slightly textured case. For TOP as well as for 
random sampling the total error is slightly lower than for the previously 
investigated CODF, e.g., for 32 grains the total error is 6'°t = 3.62% and 
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ĝt = 0.40% for random sampling and TOP, respectively. This holds for 
higher grain numbers as well. Indeed, for 1536-grain microstructures, 
randomly sampling orientation data leads to a total error of 6° = 0.53%, 
whereas using orientations generated by TOP results in an error of 
st — 0.06%. However, the relative difference between the total errors 
for the two CODFs is lower for random sampling than for TOP. 
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Figure 4.11: Total error 5' of the effective stiffness computed for microstructures with a 
unique and log-normal grain size distribution and a highly textured CODF. (a) Unique 
grain size distribution, (b) Log-normal grain size distribution. 


For the case of a log-normal GSD, compared to microstructures with 
grains of equal volume, the total error increases for both methods. For 
microstructures with 64 grains equipped with orientations randomly 
sampled from experimental data, the total error is 6'°' = 2.82% and 
tet = 3.11% for a uniform and log-normal GSD, respectively. This 
observation holds for the TOP method as well, e.g., using 64 grains leads 
to an error increase from ôt% = 0.40% for a unique GSD to ö'* = 0.45% if 
the grain sizes follow a log-normal distribution. For the TDT algorithm 
and a grain count below 768, we set the number of elementary volumes 
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for the smallest grain to eight. To account for the increased grain count, 
we increase the number of elementary volumes to twelve for 1024 and 
1536 grains in the microstructure, whereas we retain the threshold of 
5° for the misorientation computations. For our choice of parameters 
and grain numbers up to 256, we observe that the resulting error is 
close to random sampling. For instance, the total error obtained using a 
256-grain microstructure is 6'°t = 0.94% and 5' = 1.37% for orientations 
from random sampling and the TDT algorithm, respectively. When 
increasing the grain count, there seems to be a limiting accuracy that 
the TDT algorithm can reach. Indeed, the total error does not decrease 
below 1%. The source of this phenomenon needs to be investigated more 
thoroughly, and is beyond the scope of this work. 


4.3.3 Cyclic stress-strain hysteresis 


We expand our investigation into the elasto-plastic regime, focusing on 
the effect of the orientation sampling method on the cyclic stress-strain 
hysteresis of the material. As boundary condition, we use a macroscopic 
strain which follows a triangular path with an amplitude of £a = 0.7% 
and a cycle time of four seconds. To ensure a stabilized cyclic stress-strain 
hysteresis, we compute two cycles in total and use the last one as our 
quantity of interest (Schafer et al., 2019a). 

Because of the increased computational cost, we restrict the investiga- 
tions to grain counts 


G € {32, 64, 128, 256, 512, 1024} 


and use five realizations per number of grains, i.e., N = 5. We use the 
material parameters specified in Tab. 4.1. 


Uniform CODF For the uniformly distributed orientations, we asses 
the isotropy of the results. For this purpose, and a single realization, 
we compute three load cases to obtain the cyclic stress-strain hysteresis 
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in three different directions, i.e., X X-, YY- and ZZ-direction. For a 
perfectly isotropic response, the stress values would coincide for every 
considered direction. Thus, to measure the deviation from this isotropic 
result, we use the average stress values 


_ 1 
Is = z(oxxs +0yys + 022) (4.24) 


of all directions at each time step s as our reference. Then, for each 
realization n, we compute the sum of the squared relative differences 
between the stresses in every direction and the mean of all directions 
weighted by the number of stress values S, i.e., 


ghysiso _ Sn OXXs 4 a OYYs _ 1 ge 9225 _ 1 . 
|: a a as 


(4.25) 


Equation (4.25) measures the mean relative deviation in all directions 
from the mean stress value (4.24). The quantity drysiso extends the 
isotropy error defined for the stiffness, see equation (4.22), where the 
ideal isotropic case corresponds to the mean stress values in all directions. 
We compute the mean error of all realizations N = 5 by 


N 
pbys,iso — = 5 Wai (4.26) 


n=l 


to get confidence in our results. 

We use the microstructures from Sec. 4.3.2 with orientations prescribed 
by Neper, TOP and random sampling. For an increasing number of 
grains, we show the error öhysiso in Fig. 4.12a. Neper and TOP behave 
similar to each other, both lying below the error values achieved by 
random sampling. For example, the mean errors obtained from mi- 
crostructures with 32 grains are dYS*S° = 5.63%, dhYss° = 5.97% and 
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dhysiso — 1.10% for random sampling, Neper and TOP, respectively. For 
an increasing number of grains in the microstructure and randomly 
sampled orientations, the error decreases more slowly than for the 
other methods. Also, random sampling results in the highest mean 
error value of d'Ys8° = 4.49% for 1024 grains. We observe a steeper 
decrease in 5*YS/S° for an increasing number of grains for Neper and 
TOP, both lying close to each other. For example, microstructures with 
64 grains produce an error of di = 1.19% and Pysio = 1.40% 
for TOP and Neper orientations, respectively. For 1024 grains, the 
error levels are sise = 0.27% and sise = 0.16% for Neper and 


TOP orientations, respectively. 
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Figure 4.12: Hysteresis isotropy error, equation (4.25) for a varying number of grains G. 
(a) Different orientation generation methods, (b) Different considered texture coefficients 


for TOP. 


Following the procedure in Sec. 4.3.2, in addition to investigating the 
degree of isotropy, we would like to assess the ability to reproduce 
the effective mechanical response with a minimum number of grains. 
In the case of non-linear mechanical behavior, we define our ground 
truth as the stress-strain hysteresis computed for five realizations of 
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a microstructure with 10000 grains, discretized by 128° voxels. From 
these five realizations, we compute the mean stress-strain curve as 


IA 

Or,s = N 5 Or sns (4.27) 
n=1 

where or s n denotes the macroscopic stress value in direction r at a given 

time step s and N refers to the number of realizations, i.e., N = 5. For 

the considered loading directions, the resulting stress-strain hysteresis 

are shown in Fig. 4.13. We observe that the individual curves lie on top 


of each other, i.e., there is no anisotropy present. 
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Figure 4.13: Comparison of mean macroscopic stress-strain hysteresis in different loading 
directions computed using five microstructures with 10000 grains and uniformly 
distributed TOP orientations. (a) XX and YY, (b) XX and ZZ, (c) YY and ZZ. 


For each number of grains G and each realization, we compute the root 
of the mean squared relative error in each direction as 


2 
is 
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where r refers to the loading direction, i.e., r € {XX,YY, ZZ} in our 
case. Then, we compute the mean error over all considered directions 


gyet I = 2 > ohysst (4.29) 
N n=1 R r=1 = 


where R denotes the number of considered directions, i.e, R = 3 
in this case. 

Comparing öhysst for different orientation sampling methods in Fig. 4.14, 
we observe similar trends as for the isotropy error 6"Y%5°. For all meth- 
ods, the error decreases with an increasing number of grains in the 
microstructure. For a small number of grains, the error resulting from 
TOP orientations is smallest with 5°Y*8' = 1.80% and 5°Y°8t = 0.78% for 
32 and 256 grains, respectively. The error from using Neper orientations 
is higher, with dhysst — 3.93% and obysst = 0.82%. Randomly sampling 
the Haar distribution results in an error of 9.00% and 2.15% for the 
same number of grains. The error for all three methods and 1024 grains 
are 0.43%, 0.44% and 1.18% for TOP, Neper orientations and random 
sampling, respectively. 

The error ö"Y®8t for taking only the texture coefficient of rank four into 
account, is shown in Fig. 4.14b together with the previously discussed 
results for considering texture coefficients with rank four and six. We 
observe that the error when accounting solely for rank four texture 
coefficients is higher than the error produced when considering higher 
ranks. The difference is less pronounced than for 5?Ys#*°, 
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Figure 4.14: Hysteresis total error, equation (4.29) for a varying number of grains G with 
uniform CODF and unique GSD. (a) Different orientation generation methods, (b) Different 
considered texture coefficients for TOP. 


To extend our studies to a non-unique grain size distribution, we use 
microstructures with a log-normal grain size distribution. We fix the 
mean and standard deviation to mean = 1 and stdev = 0.15, respectively. 
To reduce the computational effort and because Neper and TOP provided 
the most promising results, we only consider orientations generated by 
Neper and TOP in the following. 

Fig. 4.15a shows 6°Y*#° for an increasing number of grains G. For every 
number of grains, the TOP methods provides a smaller error compared 
to Neper. For instance, using 32 grains, the error for TOP and Neper is 
dhysiso — 3.72% and sise — 12.02%, respectively. The influence of the 
underlying GSD manifests. Indeed, for both cases, the values are larger 
than for the unique grain size distribution. 

Similar to the uniform GSD, the hysteresis error closely follows the trend 
observed fo the isotropy error, see Fig. 4.15b. 
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Figure 4.15: Comparison of $"Y®8! and ö"Y®$° for varying number of grains G with uniform 
CODE and log-normal GSD. (a) Hysteresis isotropy error, equation (4.25), (b) Hysteresis 
total error, equation (4.29) 


Textured CODF For the case with mild anisotropy, we consider the 
synthetic CODF described in Sec. 4.3.2. We compute the stress-strain hys- 
teresis using five microstructures consisting of 10000 grains discretized 
by 128° voxels, see Fig. 4.4a. In accordance with the case of a uniform 
orientation distribution, we use the mean stress values of these five 
realizations as our ground truth. The resulting stress-strain hysteresis 
are shown in Fig. 4.16 for all three considered loading directions. We 
observe a slight anisotropy in Y Y-direction, whereas the stress-strain 


curves in X X- and ZZ-direction coincide. 


113 


4 Generating polycrystalline microstructures with prescribed tensorial texture coefficients 


oe XX 4+ YY =— ZZ 


@ in MPa 
5 in MPa 
5 in MPa 


| i li li li L li 1 | 1 i i li li i 1 i | L li L 
0.6 -04 -0.2 0 02 04 06 0.6 —0.4 -02 0 02 04 06 —0.6 —0.4 -02 0 02 04 06 


cin% cin% cin% 


(a) (b) (c) 


Figure 4.16: Comparison of mean macroscopic stress-strain hysteresis in different 
directions computed for a microstructure with 10 000 grains and a slightly textured CODF. 
(a) XX and YY, (b) XX and ZZ, (c) YY and ZZ. 


We show the total error to the mean stress values, i.e., equation (4.29), in 
Fig. 4.17a for the TOP method as well as for randomly sampling from 
given orientations. 

Using random orientation sampling produces a larger error for all grain 
numbers considered. Especially for a small number of grains, the TOP 
method results in a visibly smaller error than for random sampling 
the experimental data. For 32 grains, the hysteresis total error 5'Y® tot 
is 3.60% and 1.28% for random sampling and TOP, respectively. For 
random sampling, the error reduces to 0.76% for 1024 grains, which is 
close to the value achieved by TOP, with 5'Y°8t = 0.52%. 
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Figure 4.17: öhys8t for varying number of grains G with a slightly textured CODF. 
(a) Unique GSD, (b) Log-normal GSD. 


We observe similar behavior for the case of a log-normal grain size 
distribution with mean = 1 and stdev = 0.15 in Fig. 4.17b. The error 
for 32 grains increases for both kinds of orientations sampling methods, 
namely to 3.84% and 1.41% for random and TOP sampling, respectively. 
For both the unique and log-normal case, similar errors of 0.85% and 
0.44% are achieved for randomly sampling experimental orientations 
and using 1024 grains. Interestingly, the error for a log-normal GSD is 
actually smaller than for the unique GSD. For instance using the TOP 
method and 512 grains, we observe an error of 5°Y°8' = 0.66% and 
orysst = 0.49% for the unique and log-normal distributions, respectively. 
For the TDT algorithm, we observe a lower error than for random sam- 
pling when an intermediate number of grains is considered, i.e., for grain 
counts of 64, 128 and 256. For instance, the error for a microstructure 
consisting of 64 grains equipped with orientations of the TDT algorithm 
is dhysst = 2.23% and dhysst = 2.22% for 128 grains. With 64 and 128 
grains, the random sampling leads to an error of dhysst — 3.84% and 
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dhysst — 2.80%, respectively. For a larger number of grains, above 256, 
the error computed for the TDT algorithm exceeds the error of the 
randomly sampled orientations from given data. For 1024 grains, the 
error for random sampling is 5°Y°8' = 0.85% and syst = 1.00% for the 
TDT algorithm. All of the observed error values are above the errors 
obtained by TOP, e.g., using 32 grains, the error for the TOP method is 
ohysst — 1.42% whereas the TDT algorithm and random sampling lead 
to an error of 5°Y58t = 4.37% and 6'Y58t = 3.84%, respectively. 


Highly textured CODF Last but not least, we consider the GSD out- 
lined in Sec. 4.3.2 with an increased degree of anisotropy. Similar to 
the case of a slight anisotropy, we compute the stress-strain hysteresis 
for five microstructures consisting of 10 000 grains and a discretization 
of 128° voxels. As our ground truth we use the mean stress of these 
five realizations. 

For microstructures with grains having a unique grain size distribution 
and orientations from TOP or randomly sampling experimental data, 
we show the total error in Fig. 4.18a. We observe that, for all grain 
counts considered, using the TOP method results in lower error values 
compared to randomly sampling from given orientation data. For 
instance, using microstructures with 32 grains leads to a total error 
of dhys8t = 7.54% and dhysst = 17.44% for TOP and random sampling, 
respectively. Thus, we observe an increase in the total error in com- 
parison to the slightly textured CODF for both sampling methods and 
all microstructures. Indeed, for a 1024-grain microstructure equipped 
with orientations from TOP, the error increases from ®S8t = 1.42% 
for the slightly textured case to o'Y°8' = 5.24% for the case of higher 
texture. This observation holds for randomly selecting orientations from 
given orientation data, e.g., for a microstructure consisting of 256 grains 
the error for the highly textured CODF is 5*Y°8' = 10.1% whereas it is 
ohysst = 1.82% for the slightly textured case. 
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Figure 4.18: 5"Y°8' for varying number of grains G with a highly textured CODE. (a) Unique 
GSD, (b) Log-normal GSD. 


Fig. 4.18b shows the total error for the case of a log-normal GSD 
equipped with orientations from TOP, TDT and random sampling. We 
make similar observations to the slightly textured CODE  i.e., an increase 
in the total error compared to the results for the unique GSD. For instance, 
a 64-grain microstructure with TOP orientations leads to an increase in 
the total error from 6Y°8t = 5.36% for a unique GSD to o'Ys8t = 9.62% 
for a log-normal GSD. Comparing the same microstructures equipped 
with orientations from randomly sampled orientation data, the error 
increases from 8t = 13.2% to dhysst = 16.3% for a unique and a 
log-normal GSD, respectively. For lower grain counts, i.e., below 256, 
using the TDT algorithm results in similar error values as randomly 
sampling orientation data. An exception is the 32-grain microstructure, 
for which the total error value 5'Y°8' = 11.7% is close to the value 
obtained using TOP, i.e., 5*YS8' = 8.02%. For higher grain counts, i.e., 
above 256, the TDT algorithm does not decrease the total error but 
instead we observe an increase in the error values obtained, in line with 
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observations made investigating the performance for the linear elastic 
properties in Sec. 4.3.2. 


4.4 Conclusion 


In this work, we proposed to use the coefficients of a tensorial Fourier 
expansion of the crystallite orientation distribution function (Guidi et al., 
1970) to equip digital polycrystalline microstructures with crystallo- 
graphic orientations for micromechanical simulations. Our proposed 
method is based on minimizing the difference between the current 
and the prescribed tensorial Fourier, or texture, coefficients and uses a 
gradient descent scheme on the Lie group SO(3). 

We compared the proposed texture optimization for prescribing orien- 
tations (TOP) method to different state-of-the-art methods, e.g., imple- 
mented in the sophisticated microstructure generation tool Neper (Quey 
et al., 2011; 2018). In a first step, we investigated the homogenized 
stiffnesses of polycrystals for the case of a uniform and two textured 
crystallite orientation distribution functions (CODF). Subsequently, we 
extended our studies to the non-linear case, where we investigated the 
macroscopic cyclic stress-strain hysteresis of the microstructures. For 
both, the linear elastic and the non-linear plastic case, we considered 
a unique as well as a log-normal grain size distribution (GSD). By 
introducing suitable error measures we investigated and compared the 
performance of the proposed method. 

In the isotropic, linear elastic case, TOP provided better results compared 
to the Neper method and random sampling of orientations. Using TOP, 
an isotropic effective stiffness could already be achieved for small grain 
numbers. Owing to the fact that the volume fraction of each individual 
grain is explicitly accounted for, the advantage becomes more pro- 
nounced when dealing with microstructures having a log-normal grain 
size distribution. Additionally, with TOP it is possible to reproduce the 
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linear-elastic behavior of polycrystals with a unique GSD and a textured 
CODF more accurately and efficiently than via a random sampling from 
experimental orientation data. This holds as well for the case of a log- 
normal grain size distribution. Comparing to the Texture Discretization 
Technique (TDT) algorithm proposed by Melchior and Delannay (2006), 
which also considers grain size during orientation assignment, the TOP 
method performed better for both CODFs considered. Our intensive 
numerical studies revealed that the performance of the TDT algorithm, 
in our setting, critically depends on the choice of parameters, i.e., the 
misorientation value and the number of elementary volumes per grain. 
For the non-linear plastic behavior, the results of the Neper method were 
very similar to the ones provided by TOP, showing the capabilities of 
the dedicated algorithm. Although the effect was less pronounced than 
for the case of linear-elastic behavior, we observed that a underlying log- 
normal GSD results in a decreased performance for the Neper method. 
TOP, on the other hand, was able to produce similar results as for the 
unique GSD. In addition, microstructures with orientations provided 
by the TOP method allow to accurately compute the effective, non- 
linear behavior of polycrystals with an underlying texture. For both 
textured CODFs considered in this work, the defined error measures 
were lower compared to randomly sampling orientation data. The 
errors for a log-normal GSD obtained by TOP were below the ones 
obtained with random sampling or the algorithm of Melchior and 
Delannay (2006), even for a small number of grains. Investigating a 
highly textured CODF and physically non-linear visco-plastic behavior, 
we observed higher error values compared to the slightly textured 
CODF for all algorithms considered. This contrasts with the linear 
elastic case where we observed a smaller error for the highly textured 
CODF. However, using orientations generated by the TOP method 
leads to significantly smaller errors than using the other two algorithms. 
Concerning the maximum number of texture coefficients which should 
be taken into consideration, we observed that, for the linear elastic case, 
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a low tolerance and only the texture coefficient of rank four a sufficient. 
For the non-linear behavior, we observed that accounting for the texture 
coefficient of rank six is beneficial. 

As the computational effort of micromechanical studies is mainly domi- 
nated by computing the effective behavior we omitted a comparison of 
the computational performance of TOP to the other methods. 

To conclude, we showed that extracting relevant data from the CODF 
in terms of tensorial texture coefficients leads to the most flexible and 
performing method for generating crystallite orientations for digital 
representations of polycrystalline microstructures. As orientation as- 
signment is typically treated as a post-processing step in microstruc- 
ture generation, it is possible to couple the proposed algorithm with 
well-established microstructure generators (Groeber and Jackson, 2014; 
Prasad et al., 2019; Henrich et al., 2020). With this modular structure, it 
is possible to use TOP for generating polycrystalline representations for 
a variety of applications (Flipon et al., 2020b; Vajragupta et al., 2020). 
As an additional benefit, for generating orientations, the TOP method 
requires only nine (or 22 variables, depending on the highest texture 
coefficient rank considered) to be stored, contrasting with methods that 
rely on the entire experimental database. Because of this low number 
of parameters, it is possible to fuel data driven methods (Piitz et al., 
2020; Gajek et al., 2021). Additionally, as experimental data is always 
afflicted with some degree of measurement uncertainty, investigating 
the influence of the texture on the overall macroscopic response might be 
an interesting topic, i.e., via uncertainty quantification (Bandyopadhyay 
et al., 2019; Kasemer et al., 2020). 
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Chapter 5 


Identifying material parameters 
in crystal plasticity by 
Bayesian optimization ' 


5.1 Introduction 


Using suitable methods to generate representations of polycrystalline 
microstructures, see Chap. 3 and Chap. 4, we still need a compliment 
material model to predict the cyclic behavior of polycrystalline metals. 
The crystal plasticity method, outlined in Sec. 2.2.2, proves to be a 
powerful method but needs suitable model parameters to capture the 
relevant features of the cyclic deformation behavior. 

Within this chapter we propose using Bayesian optimization with 
Gaussian processes for calibrating single crystal parameters inversely 
based on polycrystalline experiments. We improve upon related 
research (Schafer et al., 2019a) by assessing the capabilities of our 
proposed framework by comparing it to a number of extremely powerful 
methods, including a genetic algorithm and modern derivative-free 
optimization schemes. This chapter is structured as follows. We 
recall Gaussian-process based Bayesian optimization in section 5.2. 


! This chapter is based on the publication by Kuhn et al. (2021) whereas minor 
typographical and formatting changes have been made for cohesion of the manuscript 
and the readers convenience 
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We investigate the performance of the proposed approach in section 
5.3, carefully studying the optimal algorithmic parameters, comparing 
to state-of-the-art genetic algorithms (Schäfer et al., 2019a) as well as 
powerful derivative-free optimization methods (Rios and Sahinidis, 2013; 
Huyer and Neumaier, 1999; 2008) and demonstrating the applicability of 
the method to large-scale polycrystalline microstructures in section 5.3.6. 


5.2 Bayesian optimization 


Bayesian optimization (BO) is an approach for solving optimization 
problems where the objective function is expensive to compute and the 
gradient of the objective function is not available. Such optimization 
problems occur, for instance, when evaluating the function corresponds 
to running a (possibly large-scale) simulation, as is common in virtual 
crash tests (Raponi et al., 2019) or for complex chemical reactions (Hase 
et al., 2018). Such problems might also be approached by automatic 
differentiation (Griewank and Walther, 2008), which might prove dif- 
ficult to integrate into an existing simulation code, or by evolutionary 
algorithms, like particle swarm algorithms (Blum and Merkle, 2008). 
However, the latter are typically limited to small spatial dimensions and 
inexpensive function evaluations. 

In contrast, Bayesian optimization constructs a surrogate model of 
the optimization problem based on probabilistic ideas, accounting for 
uncertainty by Bayesian statistics. We restrict to GAUSSIAN PROCESS 
REGRESSION as our uncertainty model. For other approaches we refer 
to Shah et al. (2014) and Kushner (1964). 

Suppose that we are concerned with the optimization problem 


f(x) — min, (5.1) 
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where A C R? is a non-empty set in d spatial dimensions whose 
membership is easily tested (for instance, a box). For (Gaussian process 
based) Bayesian optimization, we need to specify: 


1. A mean function m : A > R. 
2. A kernel (or covariance) function K : A x A > R. 


3. An ACQUISITION FUNCTION w+ u, : A + R, depending on a value 
f* € R and two functions 


pe: A+R and o:A>Rso, 


which we will refer to as (the local estimates of) the mean and stan- 
dard deviations of the function values at each x € A. 


There are specific prerequisites for sensible kernel functions. We refer 
to section 4 in Rasmussen-Williams Williams and Rasmussen (2006) for 
details, and will discuss our choice below. 

For given m, K and u, Bayesian optimization proceeds as follows. Sup- 
pose that the objective function f was evaluated at n points x1,...,£n 
in A already, i.e., the values f(a1),..., f(£n) are known. In the Bayesian 
approach, it is assumed that the vector 


fn = [f (a1), ---, f(n) E R” 


was drawn at random from a probability distribution. In Gaussian 
process regression (Mockus, 1994), the latter vector is assumed to follow 
a multivariate normal distribution with mean un € R” and ©, E€ R"*” 
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given by 
in = [m{a1),...,m(2n)]” and 
K(a1,21) K (a1, £2) K(21,2%n) 
"ON K(x2,21) K(z2,£2) ... K(xa,x,) (5.2) 
K(x2n,21) Kl&n,22) ... Kantr) 


To infer the value f(x) at some new point we assume that the joint 
distribution of the values of f over (x1,...,&n,x) is also governed by a 
multivariate normal distribution with mean un+1 and %+1 for 2n41 = 
x, as prescribed in equation (5.2). Using Bayes’ rule, the conditional 
distribution of f(x) is also Gaussian with mean 


u(x) = m(ax) + Ku(a)" En" (fn — Hn) (5.3) 
and variance 
a(x)? = K(x, x£) + Ku(2)" 2, Kunz), (5.4) 
where 
Kn: A>Rİ, x [K(a,21),K(a,22),...,K(a,an)]" , 


see section 2 in Williams and Rasmussen (2006). Provided ©, is non- 
degenerate, and as K„(x;) corresponds to the i-th row of =, for i = 
1,..., n, we immediately see that 


ul) = f(a) and ø(z;)=0 


hold, i.e., the function j interpolates the known values f(x;), and the 
corresponding standard deviation o(x;) vanishes. 
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Furthermore, u(x) serves as an estimate for the function value f(x) at 
any potentially newly sampled point x € A. 

In this work, we set the mean function identical to zero, m(x) = 0. As 
kernel function we consider the anisotropic "Matern5/2" kernel (Matern, 
2013) 


K (x, 2) = o? (1 + v58 (a,x) + > ô (a, 2?) exp (-v55 (zx, z)) (5.5) 


in terms of a positive parameter o7 and an anisotropic distance function 


D 
5 (xi, xj) = 2o mal with ij=l...nt1 66) 
d=1 


involving dilation parameters £4 for each of the dimensions d= 1... D 
of the vector x. Indeed, we cannot justify using the squared exponential 
kernel (Stein, 2012) due to the lack of sufficient smoothness of our 
objective function f. For an overview of alternative kernel functions, 
we refer to chapter 4 in Williams and Rasmussen (2006). The dilation 
parameters of the distance function 6 may be interpreted as prefactors 
used for a nondimensionalization of the variables. To determine the 
dilation parameters and o, a maximum likelihood estimation may 
be used, see chapter 6 in Bishop (2006). For each new observation 
the Gaussian process regression is updated, i.e., the parameters of the 
kernel function (5.5) are determined with respect to all previously made 
observations, continuously improving the model of f. 

We restricted to discussing Gaussian processes for our noiseless appli- 
cation of Bayesian optimization and refer to Williams and Rasmussen 
(2006) and Bishop (2006) for a more general discussion. 

As the next step, BO searches for an improved guess for the solution 
x of the optimization problem (5.1). A possible approach would be 
to ignore the estimated statistics and to minimize u, exploiting the 
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currently estimated objective values. However, this approach disregards 
the uncertainty. Indeed, it might happen that the optimum x* is located 
in regions of high uncertainty. Thus, it might be better to explore first. 
Acquisition functions provide a suitable exploration-exploitation trade- 
off by combining the currently known means and variances into a single 
function to be optimized. 

Denote by f% the lowest objective value observed so far 


fa = min f(x). 
g= 


Plugging f% and the functions u and o above into the acquisition function 
gives rise to the surrogate optimization problem 


Ups molt) — max. 
In contrast to the original optimization problem (5.1), the acquisition 
function is cheap to evaluate and gradient information is available. Also, 
the acquisition function provides a trade-off between exploration, i.e., 
decreasing the uncertainty, and exploitation, i.e., searching in the vicinity 
of sites with small objective values. Thus, the next point to investigate is 
selected by maximizing the acquisition function 


Tn+1 = ATE MAX Ufy no: 


Recall that confidence intervals of a normally distributed random vari- 
able with mean u and standard deviation o have the form 


[u— éo, u+ Eo], 


where € is a parameter which determines the probability that a measure- 
ment lies in this confidence interval. For instance, a two-sided confidence 
interval with 95% probability is obtained for € ~ 1.96. 
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For the lower confidence-bound acquisition-function (urc) (Cox and 
John, 1992), the lower bound of the confidence interval is chosen as the 
proxy for the objective function f, i.e., 


uLca(2) = -ulz) + €o(z) (5.7) 


is considered with u(x) and o(x) given by equations (5.3) and (5.4), 
respectively. This acquisition function tries to improve the currently 
known least lower bound on f. 

The quantity € may be chosen as a tuneable parameter. Indeed, for small 
&, the acquisition function tends to select points with low mean (x). 
In contrast, high values of € encourage the algorithm to choose points 
where the variance, i.e., the uncertainty, is high. We restrict to a fixed 
value for € and refer to Srinivas et al. (2010) for an adaptive choice of €. 
An alternative strategy, proposed by Mockus et al. (1978) and made 
popular by Jones et al. (1998), considers the improvement 


In: RR, f= (fài f) (5.8) 


where (-)ı = max(0,-) denotes the Macaulay bracket. If f > fă, there 
will be no improvement. For lower values of f, the improvement Z(f) 
increases linearly. Taking the expectation of the improvement (5.8) w.r.t. 
the normal distribution determined from u(x) and o(x) gives rise to the 
expected improvement ugr(z), which may be expressed analytically in 
the form 


uere) = (An(a)) te) -anaa (SE), 69) 


a(x) a(x) 


where A, (x) = fž — u(x) and ¢ as well as ® denote the standard normal 
density and distribution function with mean ju(x) and standard deviation 
a(x), see Jones et al. (1998). 
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Lizotte (2008) proposed a shift of equation (5.9) 
Anke) = fa — ulz) 8 (5.10) 


by a parameter €. This enables controlling the trade-off between ex- 
ploitation and exploration and is similar to the parameter used in uno, 
see equation (5.7). 

The discussed acquisition functions are inexpensive to evaluate, and 
gradient data is available. Thus, in addition to zeroth order optimization 
methods, for instance Monte Carlo methods (Wilson et al., 2018) or 
DIRECT (Brochu et al., 2010), gradient-based solvers like BFGS (Frazier 
and Wang, 2016; Wang et al., 2018; Picheny et al., 2016) may be used. 
Still, finding the global optimum of the acquisition function may be 
non-trivial. Indeed, even if the original function f was convex, the 
surrogate problem need not be concave. For instance, the expected im- 
provement acquisition function (5.9) has local maxima at all previously 
explored points 21,..., En- 

Our workflow for maximizing the acquisition function works as follows. 
First, we sample the acquisition function u f+ 4, at 300 points. As our 
domain of interest is a box, we rely upon the first 300 points of the 
corresponding Sobol sequence (Sobol, 1967). Then, we select those ten 
points with the highest acquisition function values, and run BFGS with 
the selected point as initial point. Finally, we select the maximum value 
obtained during those BFGS runs. 

The employed Bayesian optimization workflow is summarized in Fig. 5.1. 
Please note that we rely upon a preliminary Latin-hypercube sam- 
pling (McKay et al., 2000) to sample the first ten points within the domain 
of interest A. 
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Input: Parameter bounds, Maximum number of iterations, 
Inital observations 


Y 


| Determine kernel function parameters a 


v 
| Compute the maximum of acquisition function | 


v 
| Evaluate objective function | 


Iteration 


Maximum 
number of 
iterations 


| Output: Optimal parameter set) 


Figure 5.1: Bayesian optimization flowchart. 


5.3 Computational investigations 


5.3.1 Setup 


We wish to identify material parameters for the quenched and tempered 
high-strength steel 50CrMo4, as investigated by Schäfer et al. (2019a). 
The material features a martensitic microstructure and a hardness of 
39 HRC. To determine the crystal plasticity parameters of this material, 
macroscopic strain-driven low-cycle fatigue (LCF) experiments at differ- 
ent strain amplitudes £ were performed, see Schäfer et al. (2019a) for 
details on the experimental procedure. 

The constitutive material model described in Section 2.2.2 was imple- 
mented into a user defined material subroutine (UMAT) within the com- 
mercial finite element solver Dassault Systèmes (2018). To thoroughly 
test the optimization framework, we restrict to the simple microstructure 
shown in Fig. 5.2, consisting of 8 x 8 x 8 grains, each discretized by 2° 
elements. The different colors in Fig. 5.2 represent distinct, randomly 
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sampled, orientations. Periodic boundary conditions were used, follow- 
ing the procedure proposed by Smit et al. (1998). 


Figure 5.2: Simplified volume element used to compute the macroscopic stress strain 
hysteresis. 


The microscopic material parameters we wish to identify are the 
following: 


1. The critical resolved shear stress 7, (2.36), which we assume to be 
identical for all slip systems, i.e., T? = Te holds for all n. 

2. The parameters of the Armstrong-Frederick or the Ohno-Wang kine- 
matic hardening model A and B, see formulas (2.37) and (2.38). 


The remaining model parameters are taken from the literature according 
to Schafer et al. (2019a), see Tab. 5.1. Please note that we restrict to slip 
systems in the lath-plane, as these are primarily activated according 
to Michiuchi et al. (2009). 
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cubic stiffness Cy, = 253.1 GPa Ciz = 132.4GPa Cy, = 75.8 GPa 
flow rule Yo = 0.001s~! c= 100 

Ohno-Wang power law exponent M=8 

lattice type BCC 

slip systems {110}(111) 


Table 5.1: Parameters used for the crystal plasticity model, see Schafer et al. (2019a). 


We collect the three parameters we wish to determine into a single 
three-component vector 


x= [T:, A, B|” eR}, (5.11) 


which we would like to identify. 
Following Herrera-Solaz et al. (2014) and Schäfer et al. (2019a), we 
consider the objective function 


f'(a) = 


ae 


S 
PG — osim)?, (5.12) 


quantifying the difference between a stress-strain hysteresis of low-cycle 
fatigue (LCF) experiments and the corresponding macroscopic stress- 
strain hysteresis of a micromechanical simulation. Here, $ denotes 
the number of time steps, cy? refers to the experimentally measured 
stress in loading direction and o$™ stands for the homogenized stress in 
loading direction. 

Following the experiments, we prescribe the macroscopic loading for 
the simulations to follow a triangular waveform in y-direction, see 
Fig. 5.2. In total, we simulate two cycles, where the each cycle lasts 
four seconds, and the second cycle serves as our simulation output. We 
use time increments of 0.01 seconds for computing the homogenized 
stress response, i.e., N = 400 is used in this work. 
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Please note that we disregard strain-rate dependency of the considered 
material due to the cyclic stable behavior, see Schäfer et al. (2019a). 

In this work, we wish to minimize the mean of the error for H different 
hysteresis, i.e., we consider 


fa) = 5 PO). 6.13) 


Here, the term f”(x) corresponds to equation (5.12) and a fixed strain 
amplitude. In this article, we use the experimental data for the strain 
amplitudes £a = 0.35%, €a = 0.60% and £a = 0.90%, i.e., H = 3. The 
corresponding hysteresis at half of the lifetime are shown in Fig. 5.3. 


Gyy in MPa 
Syy in MPa 
Syy in MPa 


1 i li l! L li li 1 L li L N li li f li 1 li li li L 
—0.9 -06 —0.3 0 03 06 09 —0.9 —0.6 -03 0 03 06 09 —0.9 —0.6 -03 0 03 06 09 


ae ea, Hey ae: 
Eyy in % Eyy in % Eyy in % 


(a) (b) (c) 


Figure 5.3: Polycrystalline stress-strain hysteresis of the martensitic 50CrMo4 steel at half 
of the lifetime used for the inverse optimization problem obtained by Schäfer et al. (2019a). 
(a) €a = 0.35%, (b) €a = 0.6%, (c) ca = 0.9%. 


For a start, we study the necessary time-step size At. Indeed, exceed- 
ingly small time steps increase the computational cost, whereas too 
large time steps distort the computational results. We consider a very 
small time step of At = 0.0001 s as our reference, and investigate the 
relative deviation of the error function (5.13) for time steps of increasing 
size, see Tab. 5.2. All simulations were carried out for the Armstrong- 
Frederick kinematic hardening model with the material parameters 
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found by Schäfer et al. (2019a). Based on these results, we select the time 
step At = 0.01 s, as we consider the associated deviation of 1.26% in the 
error function negligible. 


time step At ins 0.0001 0.0005 0.001 0.005 0.01 0.05 
f in MPa (5.13) 52.97 52.93 52.88 52.57 52.30 51.53 
Relative deviation in % | — 0.08 0.17 0.76 1.26 2.70 


Table 5.2: Resulting error values f for different time steps At. 


For the parameter identification in a Bayesian optimization framework, 
diagrammatically represented in Fig. 5.1, we consider a two-stage pro- 
cess. First, we start with running BO on a large box in parameter 
space, to obtain a rough estimate for the parameters. Based on these 
results, a second optimization with tighter boundaries is carried out for 
fine-tuning the parameters. Both steps are described in further detail in 
the following subsections. 

We use the GPy implementation (GPy, 2012) as a Gaussian-process 
regressor for the method described in Section 5.2. For the acquisition 
functions and the associated optimization we used an implementation 
by the Bosch Center for Artificial Intelligence (BCAI) based on Python 
3.7, numpy (van der Walt et al., 2011) and scipy (Virtanen et al., 2020). 
We use the maximum number of function evaluations as a stopping 
criterion within the BO framework. We set this value to 160, involving 
ten random initialization points and 150 BO steps, unless otherwise 
specified. For the evolutionary algorithms, included as a benchmark, we 
use the commercial software (Dynardo, 2019). 
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5.3.2 Parameter identification with a reasonably large 
search space 


For this setting, the search space is given by the parameter bounds listed 
in Tab. 5.3. We initialize Bayesian optimization with ten points from 
a Latin-Hypercube sampling. The "Matern5/2" kernel (5.5) is chosen 
as a kernel function. The lower confidence-bound acquisition-function 
urc (5.7) is employed with an exploration margin of € = 2.0. Thus, 
over-exploitation and getting stuck in a local maximum are avoided. The 
rationale for choosing this specific acquisition function and exploration 
margin, is discussed in Section 5.3.3. 


Parameter Tin MPa AinMPa B 
Lower bound 120. 1000. 0. 
Upper bound 220. 500000. 5000. 


Table 5.3: Permissible parameter space for the large search space problem. 


First, we investigate the parameter optimization for the AF model and 
investigate the smallest error, i.e., the smallest value of f, 


fine min f(z), (5.14) 


versus function evaluations N. Fig. 5.4a shows this quantity for a single 
BO run using a maximum number of function evaluations set to 400. Ob- 
jective values above 80.0 MPa were cut off, as these are obtained during 
the first ten function evaluations in the Latin-Hypercube sampling. 
Within this single run, BO reaches a minimum value of 48.80 MPa after 
171 function evaluations. However, a comparable error of 49.9 MPa is 
reached after 81 evaluations. 

BO reaches a strictly positive objective value. This indicates that the 
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material model described in Section 2.2.2 may only approximate the 
experimentally determined stress-strain hysteresis with a non-vanishing 
error. This effect may be systematic or result from stochastic fluctuations 
certainly present for the different experimentally determined hysteresis. 
The Bayesian optimization framework relies upon randomness both 
for the initialization and the optimization of the acquisition function. 
To evaluate the influence of this randomness, we restarted Bayesian 
optimization ten times. Fig. 5.4b shows the mean of ten BO runs in terms 
of the best error encountered versus number of function evaluations. Ad- 
ditionally, we computed the 95% confidence interval (95% CI) using an 
unbiased population formula for the standard deviation and Student’s 
t-distribution (Student, 1908). The bounds of the confidence interval are 
indicated by shaded areas centered at the mean value. 

Please note that the number of evaluations used by BO is limited to 160, 
including those used for the initialization. 
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Figure 5.4: Best error versus function evaluations for the AF model using BO. (a) Single 
run, (b) Ten runs. 


BO reaches a mean error of 47.83 MPa after the allowed 160 function 
evaluations with a 95% confidence interval of [46.40 MPa, 48.74 MPa]. 
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A comparable mean error of 48.78 MPa and a 95% confidence interval 
[47.10, MPa, 48.76 MPa] is reached after 121 function evaluations, corre- 
sponding to 111 BO steps excluding evaluations used for initialization. 
Thus, we observe that repeating the BO process leads, after a sufficient 
number of iterations, to rather similar results as for a single run. 

In addition to the AF model, we investigate the Ohno-Wang model 
in this work to demonstrate the general robustness and reliability of 
Bayesian optimization in this setting. The result for one BO run with a 
maximum number of 400 evaluations is shown in Fig. 5.5a. 
Qualitatively, similar conclusions may be drawn as for the AF model. 
To address the issue of randomness we use ten BO runs with different 
initialization and evaluate the results, see Fig. 5.5b. 
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Figure 5.5: Best error versus function evaluations for the AF model using BO. (a) Single 
run, (b) Ten runs. 


The observed behavior parallels the optimization for the AF model, 
where the mean error improves significantly within the first 100 function 
evaluations and stagnates afterwards. An important difference to the AF 
model are the final mean and the confidence bounds. For all considered 
optimization runs the mean best error is 43.94 MPa, which is lower than 
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the mean reached using AF. The 95% confidence interval is computed 
as [42.71 MPa, 45.17 MPa]. The broader 95% confidence interval for the 
OW model compared to the one obtained for the AF-model optimization 
indicates a higher dispersion of BO for the OW model, i.e., a higher 
influence of the initial sampling. 

To gain insight into the shape of the energy landscape for the considered 
black-box function, we compute the difference between two consecutive 
minimum values of f* at the j-th function evaluation 


Aff=fi-ff with ij=1...N, (5.15) 


as well as the Euclidean distance between a parameter set x; (see equa- 
tion (5.11)) and the previous vector x;_,. For proper dimensioning, we 
normalize each distance parameter-wise by the parameters zi, for which 
a minimum objective value was found, see Tab. 5.4, via 


D 2 
Tid = Ti-l, 
Ill DE). 6.16) 
d 


d=1 


Here, D denotes the number of dimensions of the input vector (5.11), 
i.e., D = 3 for this work. Let us denote by 


A= 5 N AT (5.17) 
r=1 
— ie 
G= p lle -l| 6.18) 


the averaged objective values and mean distances over F total optimiza- 
tion runs at the j-th function evaluation. 
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Parameter r:inMPa A*inMPa B* 
Armstrong-Frederick 207.54 41928.27 147.83 
Ohno-Wang 173.53 77953.37 473.07 


Table 5.4: Parameters used for normalization within the large search space setting. 


First, we consider the Armstrong-Frederick kinematic hardening model, 
see Fig. 5.6a. The first ten function evaluations, corresponding to sam- 
pling of the given parameter space, allow the algorithm to find a promis- 
ing starting point. Similar to Fig. 5.4b, for the first ten steps, the mean 
change in best error (encountered so far) is very high. After these 
initial steps, BO reduces the mean parameter distance while navigating 
towards a smaller error, i.e., the already gained knowledge about f 
in the vicinity of an encountered (local) minimum is exploited and the 
uncertainty associated to this region is reduced in the process. Thereafter, 
it becomes more attractive to sample points with a higher degree of 
uncertainty, encoded by a higher variance. This exploration phase is 
indicated by the mean parameter distance oscillating around a fixed 
value at about 50 -112 function evaluations. The error is reduced fre- 
quently by small margins, up to about 112 function evaluations. After 
that, both the mean parameter distance and the mean difference in best 
error decrease further, as no further improvement of the objective value 
is made, see also Fig. 5.4b. 

After about 120 function evaluations, corresponding to 110 BO steps, the 
exploitation tendency of the considered acquisition function is apparent. 
Despite decreasing the mean parameter distance,i.e., exploiting the 
already gained knowledge about the function, the best error is not 
decreased significantly. Indeed, the used acquisition function urc» 
comes with an intrinsic trade-off between exploration and exploitation 
by sampling the point where the negative mean minus the variance 
(multiplied by €) is largest. 

For the Ohno-Wang kinematic hardening model, the considered metrics 
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are shown in Fig. 5.6b. Sampling the parameter space reduces the 
encountered best error significantly, as we have seen for the Armstrong- 
Frederick model. Compared to the AF model, shown in 5.6a, is the mean 
relative parameter distance is smaller. This indicates a higher degree of 
smoothness in the objective function. 

After the initialization phase, the error is reduced, exploiting the already 
gained knowledge of the objective function f. The mean parameter 
distance decreases, i.e., neighboring points are being evaluated. After 
50 function evaluations, the improvement of the best error gets less 
pronounced, resulting in the plateau apparent in Fig. 5.5b. After about 80 
evaluations, a step with large mean parameter distance occurs, resulting 
in a subsequent decrease of the best error. 

Interpreted in terms of the energy landscape of the objective f, BO 
investigates a neighborhood with promising points first. After exploiting 
this area (up to about 64 iterations), the exploitation-exploration trade-off 
favors reducing the uncertainty. This results in a large mean parameter 
distance at about 80 function evaluations, allowing the algorithm to find 
anew area where the objective f may be decreased further. 
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Figure 5.6: Comparison of the mean best error difference and parameter distance found 
by Bayesian optimization for the different kinematic hardening models. (a) Armstrong- 
Frederick, (b) Ohno-Wang. 


For both single optimization runs shown in Fig. 5.4a and Fig. 5.5a, the 
resulting stress-strain hysteresis is shown in Fig. 5.7 besides their exper- 
imental counterparts. The results for, both, the Armstrong-Frederick 
and Ohno-Wang model, show good agreement with the experiments 
at medium and high strain amplitudes, i.e., € = 0.6% and e = 0.9%. 
These observations may be quantified when normalizing the root of 
the mean squared distance between experiment and simulation, i.e., 
equation (5.12), by to the root of the mean squared experimental values 


f” 
N xp 
y X i= A 
For the Armstrong-Frederick model, the relative difference is 4.29% and 
6.81% for e = 0.6% and e = 0.9%, respectively. We obtain similar values 
for the Ohno-Wang model, namely 3.56% and 6.48%. For the smallest 
strain amplitude considered, i.e., € = 0.35%, the relative error is larger, 


with 22.60% and 20.41% for the Armstrong-Frederick and Ohno-Wang 
model, respectively. These errors may be reduced either by considering 


(5.19) 
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a more complex microstructure or by working with more advanced 


models. For the work at hand, we focus on the optimization method 


used for the parameter identification. 
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Figure 5.7: A comparison of stress-strain hysteresis at different loading amplitudes with 
parameters identified by BO to experimental results. (a) €a = 0.35%, (b) ca = 0.6%, 


(©) Ea 


= 0.9%. 


5.3.3 On the choice of acquisition function 


In this section, we elaborate on our choice of acquisition function and 


exploration margin. For the large search-space problem, we compare 


two different exploration margins for each of the considered acquisition 


functions ugr and urcp. For these four combinations, the mean best 


error of ten BO runs is shown in Fig. 5.8a. Except for uncg with € = 1.0, 


all combinations lead to a comparable minimum value. 
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Figure 5.8: Comparison of the mean best error using BO with different acquisition functions 
and exploration margins. (a) Entire range of evaluations, (b) Last 100 evaluations. 


Additionally, we look at the maximum best error (i.e., the worst case) 
for all ten optimization runs versus the number of function evaluations, 
shown in Fig. 5.9a. As for the behavior of the mean best error, shown 
in Fig. 5.8 for urcz with € = 1.0, the maximum best error does not 
decrease to a comparable level. The other three considered cases behave 
similarly. To examine these cases in more detail, we investigate at the 
last 100 evaluations in Fig. 5.9b. BO using uzcg with an exploration 
margin of two reaches its minimum objective value of 49.45 MPa after 
119 function evaluations for the considered case. The worst case for 
expected improvement uz; and € = 0.1, a stationary value of 51.36 
MPa, is reached after 106 function evaluations. The reached value is 
not improved using more function evaluations. For an exploration 
margin of 0.05, a stationary value of 49.31 MPa is reached for the worst 
case considered. This is slightly smaller than the error obtained using 
urcp and € = 2.0, but requires 26 additional evaluations, which is 
why we chose uzcg with € = 2.0 as our acquisition function and 


exploration margin. 
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Figure 5.9: Comparison of the worst case obtained best error using BO with different 
acquisition functions and exploration margins. (a) Entire range of evaluations, (b) Last 100 
evaluations. 


5.3.4 Fine-tuning material parameters with small 
parameter space 


In section 5.3.2, we were concerned with optimizing parameters when 
little is known about the optimum. We may work on a smaller parameter 
space if additional information is known, for instance based on expert 
knowledge or a previous optimization run for a similar problem. For 
determining the crystal plasticity parameters of the single crystal, we 
will use the bounds of Tab. 5.5, which we chose to be tighter than the 
previously considered bounds. Again, we choose our acquisition func- 
tion to be uzcpg with an exploration margin of € = 2.0, as conclusions 
very similar to those in section 5.3.3 may be drawn the fine-tuning case, 


as well. 
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Parameter Tin MPa AinMPa B 
Lower bound 120. 1000. 0. 
Upper bound 220. 200000. 1200. 


Table 5.5: Permissible parameter space for the fine-tuning problem. 


The BO workflow is initialized by ten points obtained from a Latin- 
Hypercube sampling. Then, BO is run for a fixed set of 150 function 
evaluations, summing up to a total of 160. The results of ten optimization 
runs are shown in Fig. 5.10a and 5.10b for the Armstrong-Frederick and 
Ohno-Wang model, respectively. 
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Figure 5.10: Mean best error and corresponding 95% confidence interval versus function 
evaluations of the Bayesian optimization in the fine-tuning setting. (a) Armstrong- 
Frederick, (b) Ohno-Wang. 


BO reaches a mean minimum error-value of 46.72 MPa for the AF 
model. The 95% confidence interval at end of optimization is given 
by [45.66 MPa, 47.78 MPa], which is tighter than for the results of 
the large search-space problem. After 80 evaluations, BO reaches a 
comparable mean error of 46.92 MPa with a 95% confidence interval 
of [45.76 MPa, 48.08 MPa]. 
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For the optimization of the OW kinematic hardening model, Fig. 5.10b 
shows the mean and the 95% confidence interval of the smallest error 
versus function evaluations. The mean minimum-value of 41.69 MPa 
is reached after 136 evaluations and a comparable mean error of 41.695 
MPa is reached after 73 evaluations. Furthermore, each of the optimiza- 
tion runs arrives at almost the same error, emphasized by the narrow 
95% confidence interval s, [41.6851 MPa, 41.6854 MPa] and [41.52 MPa, 
42.17 MPa] for 136 and 73 function evaluations, respectively. 

As we did for the large search-space problem, we look at, both, the mean 
difference between two consecutive best encountered values of f, and, 
the mean normed parameter distance (see equations (5.17) and (5.18)). 
As normalization, we use those parameters for which the minimum 
function value was found, see Tab. 5.6. 


Parameter TinMPa A*inMPa B* 
Armstron-Frederick 203.0 44484.21 176.0 
Ohno-Wang 186.52 63853.35 403.67 


Table 5.6: Parameters used for normalization for the fine-tuning setting. 


These metrics are shown for the AF model in Fig. 5.11a. BO system- 
atically reduces the step size, with the main change made between 
0 and 60 function evaluations. This corresponds to the observations 
made in Fig. 5.10a, where after about 60 function evaluations there is no 
substantial change in the mean best error encountered as well as for the 
95% confidence interval . Fig. 5.11a also shows the trade-off between 
exploration and exploitation. This is indicated by peaks in parameter 
distance after about 80 function evaluations without decreasing f. The 
mean distance between parameters is smaller within the more restrained 
boundaries of Tab. 5.5. This explains why a comparable or even better 
solution is achieved in shorter time than when optimizing with the 
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(larger) bounds given in Tab. 5.3. 

For the fine-tuning of the OW model, see Fig. 5.11b, similar observations 
can be made. In contrast to optimizing the AF model, the mean change 
in best error, shown in Fig. 5.10b, is almost zero after 80 function evalua- 
tions. As for the large-space problem, the mean parameter distance is 
smaller than for the case of optimizing the Armstrong-Frederick model. 
However, for the case at hand, this difference is less pronounced. 
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Figure 5.11: Comparison of mean error and parameter distance obtained with BO for 
(a) the Armstrong-Frederick and (b) the Ohno-Wang kinematic hardening models within 
the fine-tuning setting. 


5.3.5 Comparison to the state-of-the-art 


In this section, we compare the proposed Bayesian approach to alter- 
native optimization algorithms which demonstrated their power for 
similar tasks, in particular inversely identifying material parameters. In 
addition to the evolutionary approach used by Schäfer et al. (2019a), we 
investigate the Multilevel Coordinate Search (MCS) (Huyer and Neu- 
maier, 1999) and Stable Noisy Optimization by Branch and FIT (SNOB- 
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FIT) (Huyer and Neumaier, 2008). The latter two derivative-free op- 
timization methods turned out to be particularly powerful for low- 
dimensional problems in the expressive review article of Rios and Sahini- 
dis (2013). The workflow used by Schafer et al. (2019a) proceeds as 
follows. First, a Latin-Hypercube-Sampling, consisting of 200 points in 
the proposed parameter space, is performed. Subsequently, an evolution- 
ary algorithm is run on the parameter space, initialized using the most 
promising results from the sampling step. MCS and SNOBFIT do not 
require an additional initialization, and we use the recommended default 
parameters (Huyer and Neumaier, 2008; 1999) for both algorithms. 

We compare the smallest error encountered versus function evaluations 
for a single optimization run of the AF model in the large-search-space 
setting, see Fig. 5.12a, where we constrain the maximum number of 
function evaluations to 400. For a start, we observe that all algorithms 
reach a non-zero function value, in agreement to previous observations. 
The minimum error obtained by the evolutionary approach of 48.20 
MPa is reached after 332 function evaluations. A comparable value 
of 49.15 MPa is encountered after 270 evaluations. Please note that 
the sampling provides a suitable initialization, and, for the first 100 
function evaluations, there is a steady decrease in the objective value f. 
MCS reaches a similar error of 48.69 MPa after 400 evaluations, with an 
error of 49.16 MPa reached after evaluating the cost function 237 times. 
Among all considered optimization algorithms, SNOBFIT reaches the 
lowest error value of 43.72 MPa after 201 function evaluations. 

All considered algorithms, except for MCS, rely to some degree on a 
random initialization (Huyer and Neumaier, 2008; Brochu et al., 2010). 
Therefore, we wish to asses the influence of randomness on the overall 
optimization results. In Fig. 5.12b, the mean value of ten optimization 
runs is shown. Additionally, we visualize the computed 95% confidence 
interval via shaded areas, centered at the mean. For the SNOBFIT 
algorithm and these multiple runs, we set the maximum number of 
function evaluations to 160, as for the proposed Bayesian strategy. As 
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MCS is a deterministic algorithm, in all of the following we will only 
show the results of a single optimization run. 
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Figure 5.12: Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT 
and MCS for the Armstrong-Frederick kinematic hardening model in the large search 
space setting. (a) Single run, (b) Ten runs. 


After 100 evaluations, the evolutionary approach reaches a mean 
minimum function value of 56.31 MPa with a 95% confidence interval 
of [52.73 MPa, 59.90 MPa], and, after 160 evaluations arrives at a mean 
and 95% confidence interval of 53.42 MPa and [50.61 MPa, 56.23 MPa], 
respectively. The smallest mean value of 48.71 MPa is reached after 
evaluating the cost function 394 times, with a comparable value of 
49.97 MPa reached after 376 evaluations. Using the SNOBFIT algorithm 
results in a smaller mean error of 52.66 MPa at 100 function evaluations, 
but a larger 95% confidence interval of [47.14 MPa, 58.19 MPa]. The 
mean best error after 160 function evaluations is lower for SNOBFIT 
than for BO and the evolutionary algorithm, namely 47.47 MPa while 
the 95% confidence interval of [44.15 MPa, 50.78 MPa], obtained by 
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using SNOBFIT, is larger than for the other two methods. 

As our next step, we compare BO to the evolutionary approach, MCS and 
SNOBFIT for optimizing the unknown parameters of the OW kinematic 
hardening model in the large-space setting. In Fig. 5.13a, we plot the 
best error vs. iterations for each algorithm and a single optimization 
run. Bayesian optimization, the evolutionary approach and SNOBFIT 
reach similar error values of 44.21 MPa, 45.74 MPa and 43.72 MPa, 
respectively. Comparable error values are obtained by BO, SNOBFIT 
and the evolutionary framework after evaluating the cost function about 
120, 290 and 340 times, respectively. Using Multilevel Coordinate Search 
leads to a minimum best error of 48.69 MPa, which is reached after 392 
function evaluations. 

To investigate the influence of the random initialization on optimizing 
the OW model, we use the mean of ten optimization runs and visualize 
the results, with the corresponding 95% confidence interval as shaded 
areas, in Fig. 5.13b for BO, SNOBFIT and the evolutionary approach. 
Combining Latin-Hypercube sampling and the evolutionary algorithm 
produces, after 400 function evaluations, a mean error of 51.10 MPa, 
and a 95% confidence interval of [47.76 MPa, 54.40 MPa]. After 
evaluating the function 160 times, the mean best error using the 
evolutionary approach is 55.30 MPa with 95% confidence interval 
[53.19 MPa, 57.40 MPa]. 

For SNOBFIT, the mean best error of 49.77 MPa after 160 function 
evaluations lies between the results of BO and the evolutionary approach, 
whereas the 95% confidence interval of [44.73 MPa, 54.81 MPa] is larger 
than for the two other methods. 
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Figure 5.13: Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT 
and MCS for the Armstrong-Frederick kinematic hardening model in the large search 
space setting. (a) Single run, (b) Ten runs. 


In Fig. 5.14, we compare the considered algorithms for the refined 
bounds given in Tab. 5.5, for both kinematic hardening models. 

For optimizing the parameters of the AF model, SNOBFIT provides 
the lowest mean error of 44.16 MPa with [42.99 MPa, 45.34 MPa] as 
a 95% confidence interval after 160 function evaluations. We observe 
comparable values after letting SNOBFIT evaluate the function 98 times, 
e.g., 44.53 MPa for the mean error with a 95% confidence interval of 
[43.32 MPa, 45.74 MPa]. Using an evolutionary algorithm to optimize 
the AF parameters leads to a mean minimum error of 47.29 MPa, which 
is close to the mean value reached by BO. The 95% confidence interval 
at 160 function calls of the evolutionary algorithm computes to [45.47 
MPa and 49.10 MPa], which is larger than for BO or SNOBFIT. To reach 
similar mean error values of 47.50 MPa and 46.85 MPa, the evolutionary 
algorithm and BO need 153 and 87 function evaluations, respectively. 
Among the considered algorithms, MCS reaches the lowest best error, 
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already in the first 13 iterations with 50.807 MPa. This error level is 
reached within the first function evaluations. Yet, MCS falls short of 
reducing the error significantly for subsequent function evaluations. 
Indeed, after evaluating the function 160 times, the error is 48.70 MPa. 
Next, we compare the mean best error encountered when optimizing 
the parameters for the Ohno-Wang kinematic hardening model using 
the evolutionary approach, BO or SNOBFIT, see Fig. 5.14b, together with 
best error over iterations for the MCS method. Similar conclusions for 
the behavior of SNOBFIT may be drawn as for the BO approach, see 
Sec. 5.3.4. The mean best error and 95% confidence interval after 160 
function evaluations are 43.10 MPa and [42.60 MPa, 43.60 MPa], respec- 
tively. Using the evolutionary algorithm produces a mean minimum 
error of 49.31 MPa after 133 function evaluations with a 95% confidence 
interval of [48.78 MPa, 49.84 MPa]. Almost the same mean best error 
of 50.71 MPa with a 95% confidence interval of [48.79 MPa, 52.64 MPa] 
is reached after 13 function evaluations and not improving up to 109 
evaluations. For MCS, the error is, similar to optimizing the AF model, 
the lowest in the first 10 optimization steps. After evaluating the cost 
function 160 times, the error reaches a similar value as reached by BO 
with 42.29 MPa. 
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Figure 5.14: Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT 
and MCS for both kinematic hardening model in the fine-tuning setting. (a) Armstrong- 
Frederick model, (b) Ohno-Wang model. 


Let us summarize the findings of this section for each algorithm indi- 
vidually. The evolutionary algorithm consistently required the highest 
number of function evaluation among the considered algorithms to 
reach a specific error level. For the OW model and both the large and 
the small search space, the evolutionary algorithm lacked significantly 
behind the other algorithms. With the exception of the large search 
space and the AF model, MCS provided a similar error value as BO 
and the evolutionary algorithm. MCS showed inferior performance for 
the OW model and the large search space. In the case of fine tuning, 
MCS turned out to be worst for the Armstrong-Frederick model, but 
second best for Ohno-Wang. Still, there are practical benefits of the MCS 
method. First and foremost, it is a deterministic algorithm, dispensing 
with the influence of randomness. Secondly, in the fine-tuning setting 
MCS provided very good results already for the first few iterations. Thus, 
when the number of evaluations is small and prescribed beforehand, 
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MCS is the method of choice. 

SNOBFIT is good for the fine-tuning case, in particular for the AF model. 
In this case, SNOBFIT produces a similar dispersion as BO. However, 
for the large search space, the values produced by SNOBFIT are charac- 
terized by a rather large dispersion. Directly comparing SNOBFIT and 
BO, the former consistently requires more function evaluations to reach 
the same error level. 

The proposed Bayesian optimization framework produces a mean which 
is comparable to SNOBFIT in the setting of the large search space. 
However, BO is both consistently faster and characterized by a smaller 
dispersion. Only for fine-tuning AF, BO is a little worse than SNOBFIT 
in terms of the reached best error. For fine-tuning the more complex 
Ohno-Wang model, BO is fastest and leads to the best error. 

To conclude, the presented BO framework provides a strongly competi- 
tive solution, independent of the problem. Thus, it turns out to be the 
algorithm of choice for industrial applications, in particular parameter 
identification for an unknown material. 


5.3.6 Industrial-scale polycrystalline microstructure 


Previously, we determined the optimum parameters for a simplified 
volume element, see Fig. 5.2, which was selected with extensive studies 
of the proposed Bayesian optimization approach in mind. In this sec- 
tion, we investigate a more complex microstructure to demonstrate the 
capabilities of the approach. For this we use a synthetic microstructure, 
generated by the method discussed in Kuhn et al. (2020), consisting 
of 100 grains (with equal volume). Furthermore, each grain is ran- 
domly assigned a crystalline orientation, s.t. the overall orientation is 
isotropic, see Fig. 5.16a, where the coloring indicates the corresponding 
orientation. We furnish this volume element with periodic boundary 
conditions, discretize it by 64° voxels and use the FFT-based solver 
FeelMath (Fraunhofer ITWM, 2021; Wicht et al., 2020a;b) to speed up 
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computations (Rovinelli et al., 2020). Both, the Armstrong-Frederick and 
the Ohno-Wang kinematic hardening model will be considered. 
Following Sedighiani et al. (2020), we define the bounds of the feasible 
parameter space in terms of the previously found optima, with a range 
of 1.5 x x* and 0.5 x x*, see Tab. 5.7 and Tab. 5.8. We retain our choice of 
acquisition function and exploration margin. To keep the computational 
cost reasonable, we restrict the maximum number of iterations to 50. 


Parameter r.inMPa AinMPa B 
Lower bound 103. 20964. 74. 
Upper bound 310. 62892. 222. 


Table 5.7: Parameter space for the AF model following Sedighiani et al. (2020). 


Parameter r.inMPa AinMPa B 
Lower bound 93. 31927. 202. 
Upper bound 280. 95780. 806. 


Table 5.8: Parameter space for the OW model following Sedighiani et al. (2020). 


The best error versus number of evaluations is shown in Fig. 5.15a 
for the Armstrong-Fredrick model. Guided by a Latin-Hypercube- 
Sampling, the first ten function evaluations result in an error of 44.10 
MPa. This error is further reduced to 43.68 MPa after 20 additional 
function evaluations. The minimum error of 41.20 MPa is reached after 
evaluating the function 33 times. The obtained error is smaller than in 
Section 5.3.4. This may be a result of the increased complexity of the 
microstructure, which is closer to the experimental setup encompassing 
a large number of grains. Furthermore, the obtained parameter values, 
shown in Tab. 5.9 lie outside of the bounds considered in the previous 
sections, see Tab. 5.3 and Tab. 5.5. 
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Figure 5.15: Smallest error versus iterations for the optimization of kinematic hardening 
parameters and critical resolved shear stress using the microstructure with increased 
complexity. (a) Armstrong-Frederick model, (b) Ohno-Wang model. 


Similar conclusions may be drawn for the Ohno-Wang kinematic hard- 
ening model, where the smallest error versus number of function evalu- 
ations is shown in Fig 5.15b. After ten function evaluations following a 
Latin-Hypercube-Sampling, the smallest error is 62.05 MPa. After taking 
one BO step, this error is reduced to 41.35 MPa. The final minimum 
error of 38.21 MPa is reached after 27 function evaluations, i.e., after 17 
steps proposed by Bayesian optimization. 


Parameter r.inMPa AinMPa B 
Armstrong-Frederick 226.55 42758.98 116.34 
Ohno-Wang 239.50 37166.05 219.87 


Table 5.9: Resulting parameter sets for the optimization with a more complex microstruc- 
ture model. 


To ensure that the smaller error for the complex microstructure is not 
a consequence of the different parameter bounds, we return to the 
simplified microstructure, see Fig. 5.2, and run Bayesian optimization 
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on the constrained parameter spaces, see Tab. 5.7 and Tab. 5.8. For the 
Armstrong-Frederick model, BO reaches a mean error of 44.69 MPa 
an a [44.69 MPa, 44.70 MPa] 95% confidence interval after 50 function 
evaluations. After the same number of evaluations, the resulting mean 
error and 95% confidence interval are 41.78 MPa and [41.59 MPa, 41.97 
MPa], respectively. Thus, even for the constrained bounds, the error 
corresponding to the more complex microstructure shown in Fig. 5.16a 
improves upon the error obtained for the simplified structure. 
Concerning the computational cost, the simulation with a time step of 
At = 0.01s and 16 CPUs takes roughly 40 minutes for the largest strain 
amplitude of £a = 0.9%. Thus, for the chosen maximum number of 50 
function evaluations, the total Bayesian parameter identification takes 
33 hours. The minimum error of 38.21 MPa is found after a computation 
time of 18 hours in total for the Ohno-Wang kinematic hardening model. 
As for the Armstrong-Frederick model a total computation time of 22 
hours is needed to find a solution with error 41.20 MPa. 

Using a more complex synthetic microstructure, the error is further 
reduced, at the expense of an increased computational effort. Still, the 
Bayesian-optimization approach limits the required total time to less 
than a day. Using a more complex representation of the microstruc- 
ture offers also permits us to gain deeper insights into the deforma- 
tion behavior of the microstructure, for instance in terms of the accu- 
mulated plastic slip, a mesoscopic indicator for the plastic deforma- 
tion. Furthermore, the obtained optimum parameter set may enter 
further microstructure simulations. 

For £a = 0.9% and t = 8.0s, the distribution of accumulated plastic slip, 
see Wulfinghoff et al. (2013), 


Ni: (5.20) 
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is shown in Fig. 5.16b and Fig. 5.16c for the Armstrong-Frederick and 
Ohno-Wang kinematic hardening model, respectively. In both cases, the 
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Figure 5.16: Resulting distribution of the accumulated plastic + slip for €a = 0.9% at the 
end of simulation. (a) Morphology, (b) AF model, (c) OW model. 


As we used the same geometric microstructure for each of the harden- 
ing models, we observe that, independent of the employed kinematic 
hardening model, the orientation of each grain plays a crucial role in 
determining the deformation behavior. Apparently, grains with higher 
accumulated plastic slip 7 are oriented in such a way that the applied 
macroscopic load has maximum effect, thus resulting in a higher amount 
of plastic deformation. 

Moreover, the distribution of accumulated plastic slip provides some 
insight into the advantages of the OW model over the AF model in 
producing smaller errors which is not apparent from the resulting hys- 
teresis shown in Fig. 5.17. For the Ohno-Wang kinematic hardening 
model, there are more regions with high accumulated plastic slip + in 
Fig. 5.16c than for the Armstrong-Frederick model in Fig. 5.16b. This 
more pronounced plasticity seems to match the experimental behavior 
better than the results computed using the Armstrong-Frederick kine- 
matic hardening model. The relative error defined in equation (5.19) is 
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considerably lower for the small strain amplitude e = 0.35% for both 
kinematic hardening models with 11.65% and 8.65% for the Armstrong- 
Frederick and Ohno-Wang model, respectively. The relative errors for 
the medium strain amplitude are increased to 11.20% for the Armstrong- 
Frederick and 11.08% for the Ohno-Wang model, whereas the error 
for e = 0.9% is also reduced to 4.59% and 5.16% for the respective 
models. As the small strain amplitude £ = 0.35% is most interesting for 
fatigue applications, the benefits of using the more complex structure, 
see Fig. 5.16a, sinstead of the simplified microstructure, see Fig 5.2, 
become apparent. 


— Experiment — AF — OW 


900 — T T T T T T 900 — T T T T T T 900 ~ 
600 H | 600 | | 600 + 


300 H 4 a 300- J a 300b 


Syy in MPa 


li 1 li | li li li 1 li 1 | li li li 1 | li li li li 
—0.9 -06 -03 0 03 06 09 —0.9 -06 -03 0 03 06 09 —0.9 -06 -03 0 03 06 09 


Eyy in % Eyy in % Ey in % 
(a) (b) (c) 


Figure 5.17: Comparison of stress-strain hysteresis determined by BO to experimental 
results at different loading amplitudes. (a) cu = 0.35%, (b) €a = 0.6%, (c) ca = 0.9%. 


5.4 Conclusion 


This work was dedicated to identifying constitutive parameters for a 
small strain crystal plasticity constitutive law incorporating a kinematic 
hardening model to capture the behavior under cyclic loading condi- 
tions, a topic which received attention recently from, both, a scientific 
and applied perspective (Sedighiani et al., 2020; Shahmardani et al., 
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2020; Shenoy et al., 2008; Prithivirajan and Sangid, 2018; Guery et al., 
2016; Hochhalter et al., 2020; Kapoor et al., 2021; Pagan et al., 2017; 
Bandyopadhyay et al., 2020; 2019; Rovinelli et al., 2018). 

We investigated a Bayesian optimization approach for the rapid and 
flexible calibration of the single-crystal constitutive parameters entering 
micromechanical simulations. The methodology is data-driven and 
compares the predicted stress-strain hysteresis to given polycrystalline 
experimental data. More precisely, we proposed using Bayesian 
optimization with Gaussian processes to build up a surrogate model of 
the optimization problem at hand and employed the upper confidence 
bound using an exploration margin of two for selecting the next iterate. 
Our numerical investigations demonstrated that the proposed optimiza- 
tion framework reliably and efficiently determines suitable material 
parameters for our purposes. For instance, these moduli may enter a 
simulation predicting the initiation of fatigue cracks, see Schafer et al. 
(2019a). We investigated both a large and a comparatively small space 
of admissible parameters, simulating whether prior knowledge on the 
optimum is available or not. The proposed framework is able to handle 
the challenges of the large search space, and is sped up for the smaller 
admissible set. 

By investigating both the changes in parameter vector and in the 
objective value, we could gain an intuitive understanding for the 
selection strategy of Bayesian optimization, providing a suitable trade- 
off between exploitation and exploration. Furthermore, we could 
get an insight into the general energy landscape and the effects of 
using different acquisition functions. 

We compared the proposed Bayesian optimization framework to a 
representative Schäfer et al. (2019a) of the powerful class of genetic algo- 
rithms (Kapoor et al., 2021; Bandyopadhyay et al., 2020; Rovinelli et al., 
2018) as well as two derivative-free optimization schemes recommended 
in the review article by Rios and Sahinidis (2013) for low-dimensional 
problems. The computational findings of section 5.3.5 demonstrated the 
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capabilities of the Bayesian optimization framework. Compared to the 
investigated evolutionary algorithm, BO turned out to be consistently 
faster, and featured a smaller dispersion. BO outperforms MCS mainly in 
the case of the large search space, whereas, for fine tuning, BO produces a 
smaller minimum achieved error for all cases, but requires more function 
evaluations to reach a comparable error level. Compared to SNOBFIT, 
BO produces a considerably smaller dispersion when considering a large 
search-space. BO requires fewer function evaluations, and provides 
similar or better results in case of fine-tuning. 

We observed that, for fitting computed stress-strain hysteresis to their 
experimental counterparts accurately, using more complex microstruc- 
tures is necessary. In particular, the hysteresis at a low strain amplitude, 
the one most relevant to fatigue applications, using a volume element 
with 100 complex grains led to more accurate results than the simplified 
volume element with 64 cubic grains. 

To further improve the fitting quality, it appears wise to investigate 
whether increasing the complexity of the underlying material model 
or tweaking the parameters used for microstructure generation has a 
higher effects. Especially in the context of industrial applications and 
the digital twin paradigm (Tuegel et al., 2011), accurate computational 
representations of components and their microstructure are needed. 
Further research might be invested into optimizing the microstructure 
descriptors, which are used as input for synthetic microstructure 
generation (Quey and Renversade, 2018; Bourne et al., 2020; Kuhn 
et al., 2020). In particular, it may be possible to determine a set of 
minimal, thus most efficient, characteristics a microstructure has to 
possess in order to produce the desired macroscopic material behavior. 
Complementing research dedicated to microstructure representations, 
the proposed optimization scheme could be improved. Despite the 
power and flexibility we experienced in the applied context, decreasing 
the number of necessary simulation runs (via fewer function calls) 
would be very much appreciated, in particular combined with a clever 
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stopping criterion. Furthermore, research could be devoted to dedicated 
transfer strategies, i.e., building upon previous experience with simi- 
lar microstructures, material models or experimental data (Golovin 
et al., 2017). Last but not least, it appears worthwhile to extend 
the Bayesian optimization framework to uncertainty quantification, 


see Bandyopadhyay et al. (2019). 
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Summary and conclusions 


To ensure the safety of industrial components made from polycrystalline 
metals, it is imperative to assess their reliability with respect to cyclic 
loading conditions, as fatigue is still one of the major root cause of failure. 
Heterogeneities on the microscale are sources of fatigue crack initiation. 
Thus, the material’s microstructure strongly impacts the fatigue behavior 
of a component. Typically, these influences are investigated by means 
of experiments. However, fatigue experiments are time and cost inten- 
sive. Simulating the mechanical behavior on a microscopic scale allows 
predicting fatigue lifetimes, taking the microstructures influence into 
account and reducing experimental effort. 

In this work we adressed two topics required for micromechanical 
simulations of polycrystalline metals: Providing representative volume 
elements (RVEs) and identifying parameters of the crystal plasticity 
model in use. With the goal of reducing the overall runtime of the 
micromechanical simulations, we proposed two methods to create syn- 
thetic polycrystalline microstructures which reduce the required RVE 
size. We used a Bayesian optimization approach to inversely identify 
the model parameters based on polycrystalline experiments, reducing 
the required simulations by about a factor of two. 

We considered Laguerre tessellations as models for polycrystalline mi- 
crostructures in Chap. 3. Based on the formulation as an convex optimiza- 
tion problem, Bourne et al. (2020) proposed a method to quickly compute 
the weights of Laguerre tessellations, so that the computed cells have 
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prescribed volume fractions. We showed that using modern gradient- 
type solvers enables computing the weights within a few seconds. For 
all numerical examples considered, the non-monotone Barzilai-Borwein 
scheme provided the lowest run-times and therefore may serve as a 
general-purpose method. For the sake of shape regularity, it might be 
beneficial that the cells are centroidal, i.e., their seeds coincide with their 
centroids. We applied Anderson acceleration to Lloyd’s algorithm used 
by Bourne et al. (2020), which considerably speeds up the computations. 
With the algorithm at hand we were able to compute shape regular 
Laguerre tessellations which reproduced an experimentally observed 
grain size distribution. Additionally, computational cells consisting of 
multiple phases can be generated, e.g., accounting for different phase 
compositions due to a heat treatment of the metal. Extending the method 
to other features of the microstructure, e.g., sphericity of the grains, 
might be an interesting step to increase the physical realism of the 
generated synthetic microstructures. 

As a next step, in Chap. 4, we developed a method for prescribing 
crystallographic orientations to individual grains of the polycrystalline 
microstructures. We used the coefficients of a Fourier series expansion of 
the crystallite orientation distribution function, to formulate an iterative 
update procedure for the crystallographic orientations in the volume 
element. The cost function was formulated as the difference between 
prescribed and realized so-called texture coefficients. We compared 
the proposed method to other state-of-the-art approaches in terms of 
required volume element size to reproduce prescribed effective behav- 
ior. We considered, both, a uniform as well as a textured orientation 
distribution and investigated linear as well as non-linear behavior. For a 
unique grain size distribution, we observed the proposed optimization 
method to require a smaller number of grains in the volume element, 
thus, reducing the size of the RVE. Extending our investigations to the 
case of a log-normal grain size distribution, we found the proposed 
"Texture coefficient Optimization for Prescribing orientations" (TOP) 
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method shows little sensitivity to the underlying grain sizes, in contrast 
to the other considered methods. To further improve the proposed 
method, it would be interesting if faster optimization schemes, e.g., a 
Newton method, could be applied. 

With a framework, enabling the creation of high-fidelity microstructures 
at hand, we proposed a Bayesian optimization approach to determine 
suitable material model parameters in Chap. 5. More precisely, using 
an inverse optimization approach and Gaussian processes as surrogate 
models, we were able to determine single crystal plasticity parameter 
from polycrystalline experiments quickly and efficiently. Compared to the 
state-of-the-art and powerful derivative-free approaches, the Bayesian 
optimization proved to be the most robust and the fastest. Applying 
transfer-learning approaches could further improve the efficiency of 
the optimization, especially when facing new, but similar, materials for 
which the parameters are to be determined. Additionally, it might be 
interesting to investigate whether the approach can be used to determine 
single phase parameters using experimental results from a multiphase 
specimen, e.g., when it is not possible to manufacture specimens with a 
single phase. 

Taking all the described results into account, this thesis supplied novel 
methods to create microstructural representations of polycrystalline 
materials and provided a fast and reliable optimization method to param- 
eterize the underlying material model. The promoted methods in this 
thesis are expected to accelerate the development of digital twins to drive 
the digitalization of the manufacturing process of metallic components. 
In order to predict the lifetime of components under various loading 
conditions using a digital twin, suitable methods for predicting the 
mechanical behavior have to be used. The methods developed in this 
thesis enable creating micromechanical models, which might be used 
either in FE?/FE-FFT approaches (Spahn et al., 2014; Kochmann et al., 
2016) or modern, Machine Learning based methods (Gajek et al., 2021). 


165 


Bibliography 


Adams, B., Olson, T., 1998. The mesostructure - property linkage in 
polycrystals. Progress in Materials Science 43 (1), 1-87. 


Adams, B. L., Boehler, J. P., Guidi, M., Onat, E. T., 1992. Group 
theory and representation of microstructure and mechanical behavior 
of polycrystals. Journal of the Mechanics and Physics of Solids 40 (4), 
723-737. 


Allen, M. P., Tildesley, D. J., 1987. Computer Simulation of Liquids. 
Clarendon Press, Oxford. 


Alpers, A., Brieden, A., Gritzmann, P., Lyckegaard, A., Poulsen, H. F., 
2015. Generalized balanced power diagrams for 3D representations of 
polycrystals. Philosophical Magazine 95 (9), 1016-1028. 


Anderson, D. G., 1965. Iterative procedures for nonlinear integral 
equations. Journal of the ACM 12 (4), 547-560. 


Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, ]., 
Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, 
D., 1999. LAPACK Users’ Guide, 3rd Edition. Society for Industrial and 
Applied Mathematics, Philadelphia. 


Anderson, M., Grest, G., Srolovitz, D., 1989. Computer simulation of 
normal grain growth in three dimensions. Philos Magaz B 59 (3), 293-329. 


Arnaudov, N., Kolyshkin, A., Weihe, S., 2020. Micromechanical modeling 
of fatigue crack initiation in hydrogen atmosphere. Mechanics of 
Materials 149, 103557. 


167 


Bibliography 


Arts, R. J., 1993. A study of general anisotropic elasticity in rocks by 
wave propagation: Theoretical and experimental aspects. Ph.D. thesis, 
Institut francais du pétrole. 


Asaro, R. J., Needleman, A., 1985. Texture development and strain 
hardening in rate dependent polycrystals. Acta Metallurgica 33 (6), 
923-953. 


Asaro, R. J., Rice, J. R., 1977. Strain localization in ductile single crystals. 
Journal of the Mechanics and Physics of Solids 25 (5), 309-338. 


Asmussen, S., Glynn, P. W., 2007. Stochastic Simulation: Algorithms and 
Analysis. Springer, New York. 


Aurenhammer, F., 1991. Voronoi diagrams - a survey of a fundamental 
geometric data structure. ACM Computing Surveys 23, 345-405. 


Aurenhammer, F., Hoffmann, F., Aronov, B., 1998. Minkowski-Type 
Theorems and Least-Squares Clustering. Algorithmica 20, 61-76. 


Ball, J. M., Carstensen, C., 2015. Geometry of polycrystals and microstruc- 
ture. arXiv e-prints 1507.05197, 1-7. 


Bandyopadhyay, R., Prithivirajan, V., Peralta, A. D., Sangid, M. D., 2020. 
Microstructure-sensitive critical plastic strain energy density criterion 
for fatigue life prediction across various loading regimes. Proceedings of 
the Royal Society A 476 (2236), 20190766. 


Bandyopadhyay, R., Prithivirajan, V., Sangid, M. D., 2019. Uncertainty 
Quantification in the Mechanical Response of Crystal Plasticity Simula- 
tions. JOM 71 (8), 2612-2624. 


Bansal, R. K., Kubis, A., Hull, R., Fitz-Gerald, J., 2006. High-resolution 
three-dimensional reconstruction: a combined scanning electron micro- 
scope and focused ion-beam approach. Journal of Vacuum Science & 
Technology B 24 (2), 554 - 561. 


168 


Bibliography 


Barbe, F., Decker, L., Jeulin, D., Cailletaud, G., 2001. Intergranular and 
intragranular behavior of polycrystalline aggregates. Part 1: FE model. 
International Journal of Plasticity 17, 513-536. 


Bargmann, S., Klusemann, B., Markmann, J., Schnabel, J. E., Schneider, 
K., Soyarslan, C., Wilmers, J., 2018. Generation of 3D representative 
volume elements for heterogeneous materials: A review. Progress in 
Materials Science 96, 322-384. 


Bari, S., Hassan, T., 2000. Anatomy of coupled constitutive models for 
ratcheting simulation. International Journal of Plasticity 16 (3-4), 381- 
409. 


Barzilai, J., Borwein, J. M., 1988. Two-point step size gradient methods. 
IMA Journal of Numerical Analysis 8, 141-148. 


Becker, R., 1991. Analysis of texture evolution in channel die compres- 
sion. 1. Effects of grain interaction. Acta Metallurgica et Materialia 39, 
1211-1230. 


Bertram, A., 2012. Elasticity and Plasticity of Large Deformations. 
Springer, Berlin Heidelberg, Heidelberg. 


Bertram, A., Bohlke, T., Gaffke, N., Heiligers, B., Offinger, R., 2000. On 
the generation of discrete isotropic orientation distributions for linear 
elastic cubic crystals. Journal of elasticity and the physical science of 
solids 58 (3), 233-248. 


Bishop, C. M., 2006. Pattern Recognition and Machine Learning. Springer, 
New York. 


Bishop, J. F. W., 1953. VI. A theoretical examination of the plastic 
deformation of crystals by glide. The London, Edinburgh, and Dublin 
Philosophical Magazine and Journal of Science 44 (348), 51-64. 


169 


Bibliography 


Biswas, A., Prasad, M. R. G., Vajragupta, N., Kostka, A., Niendorf, T., 
Hartmaier, A., 2020. Effect of Grain Statistics on Micromechanical Mod- 
eling: The Example of Additively Manufactured Materials Examined by 
Electron Backscatter Diffraction. Advanced Engineering Materials 22 (5), 
1901416. 


Blum, C., Merkle, D., 2008. Swarm Intelligence: Introduction and 
Applications. Springer, Berlin. 


Bohlke, T., 2005. Application of the maximum entropy method in texture 
analysis. Computational Materials Science 32 (3-4), 276-283. 


Böhlke, T., 2006. Texture simulation based on tensorial Fourier coeffi- 
cients. Computers & Structures 84 (17-18), 1086-109. 


Bohlke, T., Bertram, A., 2001. Isotropic orientation distributions of cubic 
crystals. Journal of the Mechanics and Physics of Solids 49 (11), 2459- 
2470. 


Bohlke, T., Bertram, A., 2003. Crystallographic texture induced 
anisotropy in copper: An approach based on a tensorial Fourier 
expansion of the codf. In: Journal de Physique IV (Proceedings). Vol. 
105. 


Bohlke, T., Jochen, K., Kraft, O., Löhe, D., Schulze, V., 2010. Elastic 
properties of polycrystalline microcomponents. Mechanics of Materials 
42 (1), 11-23. 


Böhlke, T., Lobos, M., 2014. Representation of Hashin-Shtrikman 
bounds of cubic crystal aggregates in terms of texture coefficients with 
application in materials design. Acta Materialia 67, 324-334. 


Borbely, A., Biermann, H., Hartmann, O., 2001. FE investigation of the 
effect of particle distribution on the uniaxial stress-strain behaviour of 
particulate reinforced metal-matrix composites. Materials Science and 
Engineering: A 313 (1-2), 34-45. 


170 


Bibliography 


Bourne, D. P., Kok, P. J. J., Roper, S. M., Spanjer, W. D. T., 2020. Laguerre 
tessellations and polycrystalline microstructures: A fast algorithm for 
generating grains of given volumes. Philosophical Magazine 100 (21), 
2677-2707. 


Bourne, D. P., Roper, S. M., 2015. Centroidal Power diagrams, Lloyd’s 
algorithm and applications to optimal location problems. SIAM Journal 
on Numerical Analysis 53 (6), 2545-2569. 


Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge 
University Press, Cambridge. 


Bravais, A., 1850. Les Systemes Formes par des Points Distribues 
Regulierement sur un Plan ou Dans Lespace. Journal de l'Ecole Polytech- 
nique, XIX 1, 1-128. 


Brenner, R., 2009. Numerical computation of the response of piezoelectric 
composites using Fourier transform. Physical Review B 79 (18), 184106. 


Brieden, A., Gritzmann, P., 2012. On optimal weighted balanced 
clusterings: Gravity bodies and power diagrams. SIAM Journal on 
Discrete Mathematics 26 (2), 415-434. 


Brisard, S., Dormieux, L., 2010. FFT-based methods for the mechanics of 
composites: A general variational framework. Computational Materials 
Science 49 (3), 663-671. 


Brochu, E., Cora, V. M., De Freitas, N., 2010. A Tutorial on Bayesian 
Optimization of Expensive Cost Functions, with Application to Active 
User Modeling and Hierarchical Reinforcement Learning. arXiv preprint 
arXiv:1012.2599 . 


Broyden, C. G., 1970. The convergence of a class of double rank 
minimization algorithms: 2. The new algorithm. Journal of Mathematical 
Analysis and Applications 6, 222-231. 


171 


Bibliography 


Bunge, H.-J., 1982. Texture Analysis in Materials Science: Mathematical 
Methods. Cuvillier Verlag. 


Burdakov, O., Dai, Y.-H., Huang, N., 2019. Stabilized Barzilai-Borwein 
method. Journal of Computational Mathematics 37 (6), 916-936. 


Cailletaud, G., 1992. A micromechanical approach to inelastic behaviour 
of metals. International Journal of Plasticity 8 (1), 55-73. 


Cavallini, F., 1999. The best isotropic approximation of an anisotropic 
Hooke’s law. Bollettino di Geofisica Teorica ed Applicata 40 (1), 1-18. 


Chakraborty, A., Eisenlohr, P., 2017. Evaluation of an inverse methodol- 
ogy for estimating constitutive parameters in face-centered cubic mate- 
rials from single crystal indentations. European Journal of Mechanics- 
A/Solids 66, 114-124. 


Chen, Y., Vasiukov, D., Gélébart, L., Park, C. H., 2019. A FFT solver for 
variational phase-field modeling of brittle fracture. Computer Methods 
in Applied Mechanics and Engineering 349, 167-190. 


Chunlei, X., Nakamachi, E., Xianghuai, D., 2000. Study of texture effect 
on strain localization of BCC steel sheets. Acta Mechanica Solida Sinica 
13 (2), 95-104. 


Coleman, B. D., Noll, W., 1974. The Thermodynamics of Elastic Materials 
with Heat Conduction and Viscosity. In: The Foundations of Mechanics 
and Thermodynamics. Springer, Berlin, Heidelberg, Heidelberg, pp. 
145-156. 


Cottrell, A. H., 1954. Dislocations and Plastic Flow in Crystals. American 
journal of physics 22 (4), 242-243. 


Courtney, T. H., 2005. Mechanical behavior of materials. Waveland Press, 
Long Grove. 


172 


Bibliography 


Cox, D. D., John, S., 1992. A statistical method for global optimization. In: 
1992 IEEE International Conference on Systems, Man, and Cybernetics. 
IEEE. 


Cruzado, A., Gan, B., Jiménez, M., Barba, D., Ostolaza, K., Linaza, 
A., Molina-Aldareguia, J. M., Llorca, J., Segurado, J., 2015. Multiscale 
modeling of the mechanical behavior of IN718 superalloy based on 
micropillar compression and computational homogenization. Acta 
Materialia 98, 242-253. 


Cruzado, A., Llorca, J., Escudero, J. S., 2020. Computational Microme- 
chanics Modeling of Polycrystalline Superalloys: Application to Inconel 
718. In: Ghosh, S., Woodward, C., Przybyla, C. (Eds.), Integrated 
Computational Materials Engineering (ICME). Springer, New York, pp. 
127-163. 


Dai, Y.-H., 2012. A perfect example for the BFGS method. Mathematical 
Programming 138 (1-2), 501-530. 


Dai, Y. H., Liao, L. Z., 2002. R-linear convergence of the Barzilai and 
Borwein gradient method. IMA Journal of Numerical Analysis 22, 1-10. 


Dassault Systémes, 2018. ABAQUS. Version: 2018. 
URL https: //www.3ds.com/de/produkte-und-services/ 
simulia/produkte/abaqus/ 


Deka, D., Joseph, D. S., Ghosh, S., Mills, M. J., 2006. Crystal plasticity 
modeling of deformation and creep in polycrystalline Ti-6242. Metallur- 
gical and Materials Transactions A 37 (5), 1371-1388. 


Dembo, R. S., Eisenstat, S. C., Steihaug, T., 1982. Inexact Newton 
Methods. SIAM Journal on Numerical Analysis 19, 400-408. 


173 


Bibliography 


Dobrich, C. M., Rau, C., Krill II, C. E., 2004. Quantitative characterization 
of the three-dimensional microstructure of polycrystalline Al-Sn using 
X-ray microtomography. Metallurgical and Materials Transactions A 
35 (7), 1953 - 1961. 


Dyck, A., Böhlke, T., 2020. A micro-mechanically motivated phenomeno- 
logical yield function for cubic crystal aggregates. ZAMM-Journal 
of Applied Mathematics and Mechanics/ Zeitschrift für Angewandte 
Mathematik und Mechanik 100 (4), e202000061. 


Dynardo, 2019. optiSLang. Version: 3.3.5. 
URL https://www.dynardo.de/software/optislang.html 


Eisenlohr, P., Roters, F., 2008. Selecting a set of discrete orientations for 
accurate texture reconstruction. Computational Materials Science 42 (4), 
670-678. 


Engels, J. K., Vajragupta, N., Hartmaier, A., 2019. Parameterization of 
a non-local crystal plasticity model for tempered lath martensite using 
nanoindentation and inverse method. Frontiers in Materials 6, 247. 


Ernesti, F., Schneider, M., Bohlke, T., 2020. Fast implicit solvers for phase- 
field fracture problems on heterogeneous microstructures. Computer 
Methods in Applied Mechanics and Engineering 363, 112793. 


Eshelby, J. D., 1957. The determination of the elastic field of an ellipsoidal 
inclusion, and related problems. Proceedings of the royal society of 
London. Series A. Mathematical and physical sciences 241 (1226), 376- 
396. 


Evans, C., Pollock, S., Rebholz, L. G., Xiao, M., 2020. A proof that Ander- 
son acceleration improves the convergence rate in linearly converging 
fixed point methods (but not in those converging quadratically). SIAM 
Journal on Numerical Analysis 58 (1), 788-810. 


174 


Bibliography 


Fang, H.-R., Saad, Y., 2009. Two classes of multisecant methods for 
nonlinear acceleration. Numerical Linear Algebra with Applications 16, 
197-221. 


Farooq, H., Cailletaud, G., Forest, S., Ryckelynck, D., 2020. Crystal 
plasticity modeling of the cyclic behavior of polycrystalline aggregates 
under non-symmetric uniaxial loading: Global and local analyses. 
International Journal of Plasticity 126, 102619. 


Fedorov, F. I., 1968. Theory of Elastic Waves in Crystals. Springer Science 
& Business Media. 


Fernandez, M. L., Böhlke, T., 2019. Representation of Hashin-Shtrikman 
Bounds in Terms of Texture Coefficients for Arbitrarily Anisotropic 
Polycrystalline Materials. Journal of Elasticity 134 (1), 1-38. 


Fletcher, R., 1970. A new approach to variable metric algorithms. The 
Computer Journal 13, 317-322. 


Fletcher, R., 2005. On the Barzilai-Borwein method. In: Qi, L., Teo, K., 
Yang, X. (Eds.), Optimization and Control with Applications. Springer 
US, Boston. 


Flipon, B., Keller, C., Quey, R., Barbe, F., 2020a. A full-field crystal- 
plasticity analysis of bimodal polycrystals. International Journal of Solids 
and Structures 184, 178-192. 


Flipon, B., Keller, C., Quey, R., Barbe, F., 2020b. A full-field crystal- 
plasticity analysis of bimodal polycrystals. International Journal of Solids 
and Structures 184, 178-192. 


Francois, D., 1989. The Influence of the Microstructure on Fatigue. In: 
Moura Branco, C., Guerra Rosa, L. (Eds.), Advances in Fatigue Science 
and Technology. Springer Netherlands, Dordrecht, pp. 23-76. 


175 


Bibliography 


Fraunhofer ITWM, 2021. FeelMath. Accessed: 01. July 2021. 
URL https://www.itwm.£fraunhofer.de/de/abteilungen/ 
sms/produkte-und-leistungen/feelmath.html 


Frazier, P. I., Wang, J., 2016. Bayesian Optimization for Materials Design. 
In: Lookman, T., Alexander, F. J., Rajan, K. (Eds.), Information Science 
for Materials Discovery and Design. Springer, Cham, pp. 45-75. 


Frederick, C. O., Armstrong, P. J., 2007. A mathematical representation of 
the multiaxial Bauschinger effect. Materials at High Temperatures 24 (1), 
1-26. 


Fritzen, F., Böhlke, T., Schnack, E., 2009. Periodic three-dimensional 
mesh generation for crystalline aggregates based on Voronoi tessellations. 
Computational Mechanics 43 (5), 701-713. 


Fu, A., Zhang, J., Boyd, S., 2020. Anderson Accelerated Douglas- 
Rachford Splitting. SIAM Journal on Scientific Computing 42 (6), A3560- 
A3583. 


Fullwood, D. T., Niezgoda, S. R., Adams, B. L., Kalidindi, S. R., 2010. 
Microstructure sensitive design for performance optimization. Progress 
in Materials Science 55 (6), 477-562. 


Gajek, S., Schneider, M., Böhlke, T., 2021. An FE-DMN method for 
the multiscale analysis of short fiber reinforced plastic components. 
Computer Methods in Applied Mechanics and Engineering 384, 113952. 


Gelebart, L., Mondon-Cancel, R., 2013. Non-linear extension of FFT- 
based methods accelerated by conjugate gradients to evaluate the 
mechanical behavior of composite materials. Computational Materials 
Science 77, 430 - 439. 


Gel’fand, I. M., Minlos, R. A., Shapiro, Z. Y., 2018. Representations of 
the Rotation and Lorentz Groups and their Applications. Courier Dover 
Publications. 


176 


Bibliography 


Gianola, D. S., Eberl, C., 2009. Micro-and nanoscale tensile testing of 
materials. JOM 61 (3), 24. 


Gillner, K., Münstermann, S., 2017. Numerically predicted high cycle 
fatigue properties through representative volume elements of the 
microstructure. International Journal of Fatigue 105, 219-234. 


Göküzüm, F. S., Nguyen, L. T. K., Keip, M.-A., 2019. A multiscale FE-FFT 
framework for electro-active materials at finite strains. Computational 
mechanics 64 (1), 63-84. 


Goldfarb, D., 1970. A family of variable metric methods derived by 
variational means. Mathematics of Computation 24, 23-26. 


Golovin, D., Solnik, B., Moitra, S., Kochanski, G., Karro, J., Sculley, 
D., 2017. Google Vizier: A Service for Black-Box Optimization. In: 
Proceedings of the 23rd ACM SIGKDD International Conference on 
Knowledge Discovery and Data Mining. 


Görthofer, J., Schneider, M., Ospald, F., Hrymak, A., Böhlke, T., 2020. 
Computational homogenization of sheet molding compound composites 
based on high fidelity representative volume elements. Computational 
Materials Science 174, 109456. 


Gottstein, G., 2013. Physical foundations of materials science. Springer, 
Heidelberg. 


GPy, 2012. GPy: A Gaussian process framework in python. http: // 
github.com/SheffieldML/GPy. 


Graf, M., Kuntz, M., Autenrieth, H., Diewald, F., Miiller, R., 2021a. Phase 
field modeling and simulation of the evolution of twelve crystallographic 
martensite variants in austenitic parent grains. PAMM 20 (1), e202000121. 


Graf, M., Kuntz, M., Autenrieth, H., Diewald, F., Müller, R., 2021b. 
Simulation of martensitic microstructures in a low-alloy steel. Archive 
of Applied Mechanics 91 (4), 1641-1668. 


177 


Bibliography 


Griewank, A., Walther, A., 2008. Evaluating Derivatives: Principles 
and Techniques of Algorithmic Differentiation, 2nd Edition. SIAM, 
Philadelphia. 


Groeber, M. A., Haley, B. K., Uchic, M. D., Dimiduk, D. M., Ghosh, 
S., 2006. 3D reconstruction and characterization of polycrystalline 
microstructures using a FIB-SEM system. Materials Characterization 
57 (4), 259-273. 


Groeber, M. A., Jackson, M. A., 2014. DREAM.3D: A Digital Representa- 
tion Environment for the Analysis of Microstructure in 3D. Integrating 
Materials and Manufacturing Innovation 3 (1), 56-72. 


Gross, D., Seelig, T., 2017. Fracture Mechanics: With an Introduction to 
Micromechanics. Springer Berlin Heidelberg, Berlin. 


Guery, A., Hild, F., Latourte, F., Roux, S., 2016. Identification of crystal 
plasticity parameters using DIC measurements and weighted FEMU. 
Mechanics of Materials 100, 55-71. 


Guidi, M., Adams, B. L., Onat, E. T., 1970. Tensorial Representation of 
the Orientation Distribution Function in Cubic Polycrystals. Textures 
and Microstructures 19. 


Harder, J., 2001. FEM-simulation of the hardening behavior of FCC single 
crystals. Acta Mechanica 150 (3-4), 197-217. 


Häse, F., Roch, L. M., Kreisbeck, C., Aspuru-Guzik, A., 2018. Phoenics: A 
Bayesian optimizer for chemistry. ACS Central Science 4 (9), 1134-1145. 


Hashin, Z., Shtrikman, S., 1962a. A variational approach to the theory 
of the elastic behaviour of polycrystals. Journal of the Mechanics and 
Physics of Solids 10 (4), 343-352. 


Hashin, Z., Shtrikman, S., 1962b. On some variational principles in 
anisotropic and nonhomogeneous elasticity. Journal of the Mechanics 
and Physics of Solids 10 (4), 335-342. 


178 


Bibliography 


Hateley, J. C., Wei, H., Chen, L., 2015. Fast Methods for Computing 
Centroidal Voronoi Tessellations. Journal of Scientific Computing 63, 
185-212. 


Haupt, P., 2013. Continuum mechanics and theory of materials. Springer, 
New York. 


Helming, K., 1996. Texturapproximation durch Modellkomponenten. 
Cuvillier. 


Henrich, M., Pütz, F., Münstermann, S., 2020. A Novel Approach to 
Discrete Representative Volume Element Automation and Generation- 
DRAGen. Materials 13 (8), 1887. 


Herrera-Solaz, V., LLorca, J., Dogan, E., Karaman, I., Segurado, J., 
2014. An inverse optimization strategy to determine single crystal 
mechanical behavior from polycrystal tests: Application to AZ31 Mg 
alloy. International Journal of Plasticity 57, 1-15. 


Hershey, A. V., 1954. The Elasticity of an Isotropic Aggregate of 
Anisotropic Cubic Crystals. Journal of Applied Mechanics 21 (3), 236- 
240. 


Hielscher, R., Schaeben, H., 2008. A novel pole figure inversion method: 
specification of the MTEX algorithm. Journal of Applied Crystallography 
41 (6), 1024-1037. 


Hill, R., 1952. The elastic behaviour of a crystalline aggregate. Proceed- 
ings of the Physical Society. Section A 65 (5), 349. 


Hill, R., 1963. Elastic properties of reinforced solids: some theoretical 
principles. Journal of the Mechanics and Physics of Solids 11 (5), 357-372. 


Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. 
Journal of the Mechanics and Physics of Solids 13 (2), 89-101. 


179 


Bibliography 


Hochhalter, J., Bomarito, G., Yeratapally, S., Leser, P., Ruggles, T., 
Warner, J., Leser, W., 2020. Non-deterministic Calibration of Crystal 
Plasticity Model Parameters. In: Ghosh, S., Woodward, C., Przybyla, 
C. (Eds.), Integrated Computational Materials Engineering (ICME). 
Springer, Cham, pp. 165-198. 


Hoetzer, J., Veronika, R., Rheinheimer, W., Hoffmann, M. J., Nestler, B., 
2016. Phase-field study of pore-grain boundary interaction. J Ceram Soc 
Jpn 124 (4), 329 - 339. 


Horvath, C. D., 2021. Advanced steels for lightweight automotive 
structures. In: Mallick, P. K. (Ed.), Materials, design and manufacturing 
for lightweight vehicles. Woodhead Publishing, Sawston, pp. 39-95. 


Hull, D., Bacon, D. J., 2001. Introduction to dislocations. Butterworth- 
Heinemann, Oxford. 


Huntington, H. B., 1947. Ultrasonic measurements on single crystals. 
Physical Review 72 (4), 321. 


Hutchinson, J. W., 1976. Bounds and self-consistent estimates for creep 
of polycrystalline materials. Proceedings of the Royal Society of London. 
A. Mathematical and Physical Sciences 348 (1652), 101-127. 


Huyer, W., Neumaier, A., 1999. Global Optimization by Multilevel 
Coordinate Search. Journal of Global Optimization 14 (4), 331-355. 


Huyer, W., Neumaier, A., 2008. SNOBFIT-Stable Noisy Optimization by 
Branch and Fit. ACM Transactions on Mathematical Software (TOMS) 
35 (2), 1-25. 


Janssens, K., 2010. An introductory review of cellular automata modeling 
of moving grain boundaries in polycrystalline materials. Math Comp 
Simul 80, 1361-1381. 


180 


Bibliography 


Johnson, N. L., Kotz, S., Balakrishnan, N., 1994. 14: Lognormal 
Distributions, Continuous univariate distributions, 2nd Edition. Vol. 1 
of Wiley Series in Probability and Mathematical Statistics: Applied 
Probability and Statistics. Wiley, New York. 


Jones, D. R., Schonlau, M., Welch, W., 1998. Efficient Global Optimization 
of Expensive Black-Box Functions. Journal of Global Optimization 13 (4), 
455-492. 


Junk, M., Budday, J., Böhlke, T., 2012. On the solvability of maximum 
entropy moment problems in texture analysis. Mathematical Models 
and Methods in Applied Sciences 22 (12), 1250043. 


Kabel, M., Bohlke, T., Schneider, M., 2014. Efficient fixed point and 
Newton-Krylov solvers for FFT-based homogenization of elasticity at 
large deformations. Computational Mechanics 54 (6), 1497-1514. 


Kanit, T., Forest, S., Galliet, I., Mounoury, V., Jeulin, D., 2003. Deter- 
mination of the size of the representative volume element for random 
composites: statistical and numerical approach. International Journal of 
Solids and Structures 40 (13-14), 3647-3679. 


Kantorovich, L. V., 1948. On Newton’s method for functional equations. 
Doklady Akademii Nauk SSSR 59, 1237-1240. 


Kapoor, K., Ravi, P., Noraas, R., Park, J.-S., Venkatesh, V., Sangid, 
M. D., 2021. Modeling Ti-6Al-4V using crystal plasticity, calibrated 
with multi-scale experiments, to understand the effect of the orientation 
and morphology of the a and £ phases on time dependent cyclic loading. 
Journal of the Mechanics and Physics of Solids 146, 104192. 


Kasemer, M., Falkinger, G., Roters, F., 2020. A numerical study of 
the influence of crystal plasticity modeling parameters on the plastic 
anisotropy of rolled aluminum sheet. Modelling and Simulation in 
Materials Science and Engineering 28 (8), 085005. 


181 


Bibliography 


Kawasaki, K., Nagai, T., Nakashima, K., 1989. Vertex models for two- 
dimensional grain growth. Philos Magaz B 60 (3), 399 - 421. 


Kehrer, L., Wood, J. T., Böhlke, T., 2020. Mean-field homogenization of 
thermoelastic material properties of a long fiber-reinforced thermoset 
and experimental investigation. Journal of Composite Materials 54 (25), 
3777-3799. 


Kitagawa, J., Mérigot, Q., Thibert, B., 2019. Convergence of a Newton 
algorithm for semi-discrete optimal transport. Journal of the European 
Mathematical Society 21 (9), 2603-2651. 


Kochmann, J., Wulfinghoff, S., Reese, S., Mianroodi, J. R., Svendsen, B., 
2016. Two-scale FE-FFI-and phase-field-based computational modeling 
of bulk microstructural evolution and macroscopic material behavior. 
Computer Methods in Applied Mechanics and Engineering 305, 89-110. 


Kocks, U. F., Tomé, C. N., Wenk, H.-R., 2000. Texture and Anisotropy: 
Preferred Orientations in Polycrystals and their Effect on Materials 
Properties. Cambridge University Press, Cambridge. 


Korte, S., Ritter, J., Jiao, C., Midgley, P. A., Clegg, W. J., 2011. Three- 
dimensional electron backscattered diffraction analysis of deformation 
in MgO micropillars. Acta Materialia 59 (19), 7241-7254. 


Kouznetsova, V., Brekelmans, W. A. M., Baaijens, F. P. T., 2001. 
An approach to micro-macro modeling of heterogeneous materials. 
Computational Mechanics 27, 37-48. 


Krawietz, A., 1999. Parallel versus Conventional Elastoplasticity. Tech- 
nische Mechanik 19 (4), 279-288. 


Krill III, C., Chen, L.-Q., 2002. Computer simulation of 3-D grain growth 
using a phase-field model. Acta Mater 50 (12), 3059 — 3075. 


182 


Bibliography 


Kroner, E., 1958. Berechnung der elastischen Konstanten des Vielkristalls 
aus den Konstanten des Einkristalls. Zeitschrift für Physik 151 (4), 504- 
518. 


Krupp, U., 2007. Fatigue crack propagation in metals and alloys: 
microstructural aspects and modelling concepts. John Wiley & Sons. 


Kubis, A. J., Shiflet, G. J., Hull, R., Dunn, D. N., 2004. Focused ion- 
beam tomography. Metallurgical and Materials Transactions A 37 (7), 
1935-1943. 


Kühlbach, M., Gottstein, G., Barrales-Mora, L., 2016. A statistical 
ensemble cellular automaton microstructure model for primary recrys- 
tallization. Acta Mater 107, 366 — 376. 


Kuhn, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 2020. Fast 
methods for computing centroidal Laguerre tessellations for prescribed 
volume fractions with applications to microstructure generation of 
polycrystalline materials. Computer Methods in Applied Mechanics 
and Engineering 369, 113175. 


Kuhn, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 2022. Gener- 
ating polycrystalline microstructures with prescribed tensorial texture 
coefficients. Computational Mechanics 70 (3), 639-659. 

Kuhn, J., Spitz, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 
2021. Identifying material parameters in crystal plasticity by Bayesian 
optimization. Optimization and Engineering , 1-35. 


Kushner, H. J., 1964. A New Method of Locating the Maximum Point of 
an Arbitrary Multipeak Curve in the Presence of Noise. Journal of Basic 
Engineering 86, 97-106. 


Lautensack, C., 2008. Fitting three-dimensional Laguerre tessellations to 
foam structures. Journal of Applied Statistics 35 (9), 985-995. 


183 


Bibliography 


Lautensack, C., Zuyev, S., 2008. Random Laguerre tesselations. Advances 
in Applied Probability 40, 630-650. 


Lebensohn, R., Rollett, A. D., 2020. Spectral methods for full-field 
micromechanical modelling of polycrystalline materials. Computational 
Materials Science 173, 109336. 


Lebensohn, R. A., Tomé, C. N., neda, P. P. C., 2007. Self-consistent 
modelling of the mechanical behaviour of viscoplastic polycrystals 
incorporating intragranular field fluctuations. Philosophical Magazine 
87 (28), 4287-4322. 


Lehto, P., Romanoff, J., Remes, H., Sarikka, T., 2016. Characterisation 
of local grain size variation of welded structural steel. Welding in the 
World 60 (4), 673 - 688. 


Lévy, B., 2015. A numerical algorithm for L? semi-discrete optimal 
transport in 3D. ESAIM: Mathematical Modelling and Numerical 
Analysis 49 (6), 1693-1715. 


Liu, W., Lian, J., Aravas, N., Miinstermann, S., 2020a. A strategy for 
synthetic microstructure generation and crystal plasticity parameter cali- 
bration of fine-grain-structured dual-phase steel. International Journal 
of Plasticity 126, 102614. 


Liu, W., Lian, J., Aravas, N., Miinstermann, S., 2020b. A strategy for 
synthetic microstructure generation and crystal plasticity parameter cali- 
bration of fine-grain-structured dual-phase steel. International Journal 
of Plasticity 126, 102614. 


Liu, W. B., Dai, Y. H., 2001. Minimization algorithms based on supervisor 
and searcher cooperation: I - faster and robust gradient algorithms for 
minimization problems with stronger noises. Journal of Optimization 
Theory and Applications 111, 359-379. 


184 


Bibliography 


Liu, Y., Wang, W., Lévy, B., Sun, F, Yan, D.-M., Lu, L., Yang, C., 2016. 
On Centroidal Voronoi Tessellation — Energy Smoothness and Fast 
Computation. ACM Transaction on Graphics (TOG) 28 (4), article 244. 


Lizotte, D. J., 2008. Practical Bayesian Optimization. Ph.D. thesis, 
University of Alberta. 


Lobos, M., Bohlke, T., 2015. Materials design for the anisotropic linear 
elastic properties of textured cubic crystal aggregates using zeroth-, 
first-and second-order bounds. International Journal of Mechanics and 
Materials in Design 11 (1), 59-78. 


Luther, T., Könke, C., 2009. Polycrystal models for the analysis of 
intergranular crack growth in metallic materials. Engineering Fracture 
Mechanics 76 (15), 2332-2343. 


Lyckegaard, A., Lauridsen, E. M., Ludwig, W., Fonda, R. W., Poulsen, 
H. F., 2011. On the use of Laguerre tessellations for representations of 
3D grain structures. Advanced Engineering Materials 13 (3), 165-170. 


Mäkinen, J., 2008. Rotation manifold SO(3) and its tangential vectors. 
Computational Mechanics 42 (6), 907-919. 


Malitsky, Y., Mishchenko, K., 2019. Adaptive gradient descent without 
descent. arXiv e-prints 1910.09529, 1-20. 


Matern, B., 2013. Spatial Variation. Vol. 36. Springer, New York. 


Matous, K., Geers, M. G. D., Kouznetsova, V. G., Gillman, A., 2017. 
A review of predictive nonlinear theories for multiscale modeling of 
heterogeneous materials. Journal of Computational Physics 330, 192-220. 


McKay, M. D., Beckman, R. J., Conover, W. J., 2000. A Comparison of 
Three Methods for Selecting Values of Input Variables in the Analysis of 
Output From a Computer Code. Technometrics 42 (1), 55-61. 


185 


Bibliography 


Melchior, M. A., Delannay, L., 2006. A texture discretization technique 
adapted to polycrystalline aggregates with non-uniform grain size. 
Computational Materials Science 37 (4), 557-564. 


Mérigot, Q., 2011. A Multiscale Approach to Optimal Transport. 
Computer Graphics Forum 30, 1583-1592. 


Michiuchi, M., Nambu, S., Ishimoto, Y., Inoue, J., Koseki, T., 2009. 
Relationship between local deformation behavior and crystallographic 
features of as-quenched lath martensite during uniaxial tensile deforma- 
tion. Acta Materialia 57 (18), 5283-5291. 


Miehe, C., Schotte, J., Schröder, J., 1999. Computational micro-macro 
transitions and overall moduli in the analysis of polycrystals at large 
strains. Computational Materials Science 16 (1-4), 372-382. 


Mockus, J., 1994. Application of Bayesian approach to numerical 
methods of global and stochastic optimization. Journal of Global 
Optimization 4 (4), 347-365. 


Mockus, J., Tiesis, V., Zilinskas, A., 1978. The application of Bayesian 
methods for seeking the extremum. In: Dixon, L. C. W., Szegö, G. P. 
(Eds.), Towards Global Optimisation 2. North Holland, Amsterdam, pp. 
117-129. 


Molinari, A., Canova, G. R., Ahzi, S., 1987. A self consistent approach 
of the large deformation polycrystal viscoplasticity. Acta Metallurgica 
35 (12), 2983-2994. 


Morawiec, A., 2003. Orientations and Rotations. Springer. 


Moulinec, H., Suquet, P., 1994. A fast numerical method for computing 
the linear and nonlinear mechanical properties of composites. Comptes 
rendus de l’Académie des sciences. Série II. Mécanique, physique, 
chimie, astronomie. . 


186 


Bibliography 


Moulinec, H., Suquet, P., 1998. A numerical method for computing the 
overall response of nonlinear composites with complex microstructure. 
Computer methods in applied mechanics and engineering 157 (1-2), 
69-94. 


MTex, 2017. MTex Documentation. https: / /mtex- 
toolbox.github.io/ODFExport.html. 


Miiller, V., Kabel, M., Andra, H., Bohlke, T., 2015. Homogenization 
of linear elastic properties of short-fiber reinforced composites—A 
comparison of mean field and voxel-based methods. International 
Journal of Solids and Structures 67, 56-70. 


Nabarro, F. R. N., 1947. Dislocations in a simple cubic lattice. Proceedings 
of the Physical Society (1926-1948) 59 (2), 256. 


Nakamura, T., Suresh, S., 1993. Effects of thermal residual stresses 
and fiber packing on deformation of metal-matrix composites. Acta 
Metallurgica et Materialia 41 (6), 1665-1681. 


Natkowski, E., , Sonnweber-Ribic, P., Miinstermann, S., 2021a. Determi- 
nation of fatigue lifetimes with a micromechanical short crack model 
for the high-strength steel SAE 4150. International Journal of Fatigue , 
106621. 


Natkowski, E., Durmaz, A. R., Sonnweber-Ribic, P., Miinstermann, S., 
2021b. Fatigue lifetime prediction with a validated micromechanical 
short crack model for the ferritic steel EN 1.4003. International Journal 
of Fatigue , 106418. 


Nemat-Nasser, S., Hori, M., 1993. Micromechanics: overall properties of 
heterogeneous materials. Elsevier, Amsterdam. 


Nesterov, Y., 2004. Introductory Lectures on Convex Optimization: 
A Basic Course. Mathematics and its applications. Kluwer Academic 
Publishers, Springer US. 


187 


Bibliography 


Nestler, B., 2005. A 3D parallel simulator for crystal growth and 
solidification in complex alloy systems. J Cryst Growth 275 (1), 273-278. 


Nocedal, J., 1980. Updating Quasi-Newton Matrices With Limited 
Storage. Mathematics of Computation 35 (151), 773-782. 


Nocedal, J., Wright, S. J., 1999. Numerical Optimization. Springer, New 
York. 


Nolze, G., Hielscher, R., 2016. Orientations—perfectly colored. Journal of 
Applied Crystallography 49 (5), 1786-1802. 


Ohno, N., Wang, J.-D., 1993. Kinematic hardening rules with critical 
state of dynamic recovery, part I: formulation and basic features for 
ratchetting behavior. International Journal of Plasticity 9 (3), 375-390. 


Orowan, E., 1934. Zur kristallplastizität. i. Zeitschrift für Physik 89 (9-10), 
605-613. 


Pagan, D. C., Shade, P. A., Barton, N. R., Park, J.-S., Kenesei, P., Menasche, 
D. B., Bernier, J. V., 2017. Modeling slip system strength evolution in 
Ti-7Al informed by in-situ grain stress measurements. Acta Materialia 
128, 406-417. 


Peierls, R., 1940. The size of a dislocation. Proceedings of the Physical 
Society (1926-1948) 52 (1), 34. 


Peirce, D., Asaro, R. J., Needleman, A., 1983. Material rate dependence 
and localized deformation in crystalline solids. Acta Metallurgica 31 (12), 
1951-1976. 


Petrich, L., Stanék, J., Wang, M., Westhoff, D., Heller, L., Sittner, P, 
Krill IIL, C. E., Beneš, V., Schmidt, V., 2019. Reconstruction of Grains 
in Polycrystalline Materials From Incomplete Data Using Laguerre 
Tessellations. Microscopy and Microanalysis 25, 743-752. 


188 


Bibliography 


Picheny, V., Gramacy, R. B., Wild, S., Le Digabel, S., 2016. Bayesian 
optimization under mixed constraints with a slack-variable augmented 
Lagrangian. In: Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., Garnett, 
R. (Eds.), Advances in Neural Information Processing Systems. Vol. 29. 
Curran Associates, Inc. 


Polanyi, M., 1934. Uber eine Art Gitterstörung, die einen Kristall 
plastisch machen könnte. Zeitschrift für Physik 89 (9), 660-664. 


Polyak, B. T., 1987. Introduction to Optimization. Optimization Software, 
Inc., New York. 


Prager, W., 1949. Recent Developments in the Mathematical Theory of 
Plasticity. Journal of Applied Physics 20 (3), 235-241. 


Prasad, M. R. G., Vajragupta, N., Hartmaier, A., 2019. Kanapy: A Python 
package for generating complex synthetic polycrystalline microstruc- 
tures. Journal of Open Source Software 4 (43), 1732. 


Prithivirajan, V., Sangid, M. D., 2018. The role of defects and critical pore 
size analysis in the fatigue response of additively manufactured in718 
via crystal plasticity. Materials & Design 150, 139-153. 


Pütz, F., Henrich, M., Fehlemann, N., Roth, A., Münstermann, S., 2020. 
Generating Input Data for Microstructure Modelling: A Deep Learning 
Approach Using Generative Adversarial Networks. Materials 13 (19), 
4236. 


Quey, R., Dawson, P., Barbe, F., 2011. Large-scale 3D random polycrystals 
for the finite element method: Generation, meshing and remeshing. 
Computer Methods in Applied Mechanics and Engineering 200 (17-20), 
1729-1745. 


189 


Bibliography 


Quey, R., Renversade, L., 2018. Optimal polyhedral description of 3D 
polycrystals: Method and application to statistical and synchrotron 
X-ray diffraction data. Computer Methods in Applied Mechanics and 
Engineering 330, 308-333. 


Quey, R., Villani, A., Maurice, C., 2018. Nearly uniform sampling of 
crystal orientations. Journal of Applied Crystallography 51 (4), 1162- 
1173. 


Radakrishnan, B., Zacharia, T., 1995. Monte Carlo simulation of grain 
boundary pinning in the weld heat-affected zone. Metall Mater Trans A 
26 (8), 2123-2130. 


Raponi, E., Bujny, M., Olhofer, M., Aulig, N., Boria, S., Duddeck, F., 2019. 
Kriging-assisted topology optimization of crash structures. Computer 
Methods in Applied Mechanics and Engineering 348, 730-752. 


Reuss, A., 1929. Berechnung der Fließgrenze von Mischkristallen 
auf Grund der Plastizitätsbedingung für Einkristalle. ZAMM-Journal 
of Applied Mathematics and Mechanics/ Zeitschrift für Angewandte 
Mathematik und Mechanik 9 (1), 49-58. 


Richard, H. A., Bürgel, R., Riemer, A., 2014. Zyklische Belastung. 
In: Bürgel, R., Richard, H. A., Riemer, A. (Eds.), Werkstoffmechanik. 
Springer Vieweg, Wiesbaden, pp. 92-115. 


Rieger, F., Böhlke, T., 2015. Microstructure based prediction and 
homogenization of the strain hardening behavior of dual-phase steel. 
Archive of Applied Mechanics 85 (9), 1439-1458. 


Rios, L. M., Sahinidis, N. V., 2013. Derivative-free optimization: a review 
of algorithms and comparison of software implementations. Journal of 
Global Optimization 56 (3), 1247-1293. 


190 


Bibliography 


Roe, R.-J., 1965. Description of crystallite orientation in polycrystalline 
materials. III. General solution to pole figure inversion. Journal of 
Applied Physics 36 (6), 2024-2031. 


Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D., Bieler, T., Raabe, D., 
2010. Overview of constitutive laws, kinematics, homogenization and 
multiscale methods in crystal plasticity finite-element modeling: theory, 
experiments, applications. Acta Materialia 58, 1152-1211. 


Rovinelli, A., Proudhon, H., Lebensohn, R. A., Sangid, M. D., 2020. 
Assessing the reliability of fast Fourier transform-based crystal plasticity 
simulations of a polycrystalline material near a crack tip. International 
Journal of Solids and Structures 184, 153-166. 


Rovinelli, A., Sangid, M. D., Proudhon, H., Guilhem, Y., Lebensohn, 
R. A., Ludwig, W., 2018. Predicting the 3D fatigue crack growth rate 
of small cracks using multimodal data via Bayesian networks: In-situ 
experiments and crystal plasticity simulations. Journal of the Mechanics 
and Physics of Solids 115, 208-229. 


Rycroft, C. H., 2009. Voro++: A three-dimensional Voronoi cell library in 
C++. Chaos 19, 041111. 


Saad, Y., Schultz, M. H., 1986. GMRES: A generalized minimal residual 
algorithm for solving nonsymmetric linear systems. SIAM Journal on 
Scientific Computing 7, 856-869. 


Saluja, R., Narayanan, R., Das, S., 2012. Cellular automata finite element 
(CAFE) model to predict the forming of friction stir welded blanks. 
Comput Mat Sci 58, 87-100. 


Schäfer, B. J., Song, X., Sonnweber-Ribic, P., ul Hassan, H., Hartmaier, A., 
2019a. Micromechanical Modelling of the Cyclic Deformation Behavior 
of Martensitic SAE - A Comparison of Different Kinematic Hardening 
Models. Metals 9 (3), 368. 


191 


Bibliography 


Schäfer, B. J., Sonnweber-Ribic, P., ul Hassan, H., Hartmaier, A., 2019b. 
Micromechanical modeling of fatigue crack nucleation around non- 
metallic inclusions in martensitic high-strength steels. Metals 9 (12), 
1258. 


Schemmann, M., Gorthofer, J., Seelig, T., Bohlke, A. H. T., 2018. 
Anisotropic meanfield modeling of debonding and matrix damage in 
SMC composites. Composites Science and Technology 161, 143-158. 


Schneider, M., 2019. On the Barzilai-Borwein basic scheme in FFT-based 
computational homogenization. International Journal for Numerical 
Methods in Engineering 118, 482-494. 


Schneider, M., 2021. A review of nonlinear FFT-based computational 
homogenization methods. Acta Mechanica 232, 2051-2100. 


Schneider, M., Josien, M., Otto, F., 2021. Representative volume elements 
for matrix-inclusion composites—a computational study on periodizing 
the ensemble. arXiv preprint arXiv:2103.07627 . 


Sedighiani, K., Diehl, M., Traka, K., Roters, F., Sietsma, J., Raabe, D., 2020. 
An efficient and robust approach to determine material parameters of 
crystal plasticity constitutive laws from macro-scale stress-strain curves. 
International Journal of Plasticity , 102779. 


Shah, A., Wilson, A., Ghahramani, Z., 2014. Student-t Processes as 
Alternatives to Gaussian Processes. In: Kaski, S., Corander, J. (Eds.), 
Proceedings of the Seventeenth International Conference on Artificial 
Intelligence and Statistics. Proceedings of Machine Learning Research. 


Shahmardani, M., Vajragupta, N., Hartmaier, A., 2020. Robust Op- 
timization Scheme for Inverse Method for Crystal Plasticity Model 
Parametrization. Materials 13 (3), 735. 


Shanno, D. F., 1970. Conditioning of quasi-Newton methods for function 
minimization. Mathematics of Computation 24, 647-650. 


192 


Bibliography 


Shantraj, P., Eisenlohr, P., Diehl, M., Roters, F., 2015. Numerically robust 
spectral methods for crystal plasticity simulations of heterogeneous 
materials. International Journal of Plasticity 66, 31-45. 


Shenoy, M., Tjiptowidjojo, Y., McDowell, D., 2008. Microstructure- 
sensitive modeling of polycrystalline IN 100. International Journal of 
Plasticity 24 (10), 1694-1730. 


Simo, J. C., Hughes, T. J. R., 2006. Computational Inelasticity. Vol. 7. 
Springer, New York. 


Smit, R. J. M., Brekelmans, W. A. M., Meijer, H. E. H., 1998. Prediction of 
the mechanical behavior of nonlinear heterogeneous systems by multi- 
level finite element modeling. Computer Methods in Applied Mechanics 
and Engineering 155 (1-2), 181-192. 


Sobol, I. M., 1967. On the point distribution in a cube and their 
approximate evaluation of integrals. USSR Computational Mathematics 
and Mathematical Physics 7 (4), 86-112. 


Spahn, J., Andrä, H., Kabel, M., Müller, R., 2014. A multiscale approach 
for modeling progressive damage of composite materials using fast 
Fourier transforms. Computer Methods in Applied Mechanics and 
Engineering 268, 871-883. 


Spettl, A., Wertz, T., Krill II, C. E., Schmidt, V., 2014. Parametric Rep- 
resentation of 3D Grain Ensembles in Polycrystalline Microstructures. 
Journal of Statistical Physics 154, 913-928. 


Spowart, J. E., Mullens, H. E., Puchalla, B. T., 2003. Collecting and 
analyzing microstructures in three dimensions: a fully automated 
approach. JOM 55 (10), 35-37. 


193 


Bibliography 


Srinivas, N., Krause, A., Kakade, S. M., Seeger, M., 2010. Gaussian 
Process Optimization in the Bandit Setting: No Regret and Experimental 
Design. In: Fürnkranz, J., Joachims, T. (Eds.), Proceedings of the 
27th International Conference on International Conference on Machine 
Learning. Omnipress, Madison. 


Stein, M. L., 2012. Interpolation of Spatial Data: Some Theory for Kriging. 
Springer, New York. 


Stewart, G. W., 1980. The Efficient Generation of Random Orthogonal 
Matrices with an Application to Condition Estimators. SIAM Journal on 
Numerical Analysis 17 (3), 403-409. 


Student, 1908. The probable error of a mean. Biometrika , 1-25. 


Tallman, A. E., Swiler, L. P., Wang, Y., McDowell, D. L., 2020. Uncer- 
tainty propagation in reduced order models based on crystal plasticity. 
Computer Methods in Applied Mechanics and Engineering 365, 113009. 


Taylor, C. J., Kriegman, D. J., 1994. Minimization on the Lie group SO(3) 
and related manifolds. Yale University 16 (155), 6. 


Taylor, G. I., 1934. The mechanism of plastic deformation of crystals. Part 
I. - Theoretical. Proceedings of the Royal Society of London. Series A, 
Containing Papers of a Mathematical and Physical Character 145 (855), 
362-387. 


Teferra, K., Rowenhorst, D. J., 2018. Direct parameter estimation for 
generalised balanced power diagrams. Philosophical Magazine Letters 
98 (2), 79-87. 


Torquato, S., 2002. Random Heterogeneous Materials: Microstructure 
and Macroscopic Properties. Springer, New York. 


Toth, A., Kelley, C., 2015. Convergence analysis for Anderson accelera- 
tion. SIAM Journal on Numerical Analysis 53 (2), 805-819. 


194 


Bibliography 


Toth, L. S., Van Houtte, P., 1970. Discretization Techniques for Orienta- 
tion Distribution Functions. Textures and Microstructures 19. 


Tran, A., Wildey, T., 2021. Solving Stochastic Inverse Problems for 
Property-Structure Linkages Using Data-Consistent Inversion and 
Machine Learning. JOM 73 (1), 72-89. 


Truesdell, C., Noll, W., 2004. The non-linear field theories of mechanics. 
In: Antman, S. S. (Ed.), The Non-Linear Field Theories of Mechanics. 
Springer Berlin, Heidelberg, Heidelberg, pp. 1-579. 


Tu, X., Shahba, A., Shen, J., Ghosh, S., 2019. Microstructure and prop- 
erty based statistically equivalent RVEs for polycrystalline-polyphase 
aluminum alloys. International Journal of Plasticity 115, 268-292. 


Tuegel, E. J., Ingraffea, A. R., Eason, T. G., Spottswood, S. M., 2011. 
Reengineering Aircraft Structural Life Prediction Using a Digital Twin. 
International Journal of Aerospace Engineering 2011. 


Vajragupta, N., Maassen, S., Clausmeyer, T., Brands, D., Schroder, 
J., Hartmaier, A., 2020. Micromechanical Modeling of DP600 steel: 
From Microstructure to The Sheet Metal Forming Process. Procedia 
Manufacturing 47, 1540-1547. 


Vajragupta, N., Wechsuwanmanee, P., Lian, J., Sharaf, M., Miinstermann, 
S., Ma, A., Hartmaier, A., Bleck, W., 2014. The modeling scheme 
to evaluate the influence of microstructure features on microcrack 
formation of DP-steel: The artificial microstructure model and its 
application to predict the strain hardening behavior. Computational 
materials science 94, 198-213. 


van der Walt, S., Colbert, S. C., Varoquaux, G., 2011. The NumPy array: 
a structure for efficient numerical computation. Computing in Science & 
Engineering 13 (2), 22-30. 


195 


Bibliography 


Virtanen, P., Gommers, R., et al., T. E. O., 2020. SciPy 1.0: fundamental 
algorithms for scientific computing in Python. Nature methods 17 (3), 
261-272. 


Voigt, W., 1889. Ueber die Beziehung zwischen den beiden Elastic- 
itätsconstanten isotroper Körper. Annalen der physik 274 (12), 573-587. 


Silhavy, M., 1997. The Mechanics and Thermodynamics of Continuous 
Media. Springer, Heidelberg. 


Wang, Z., Gehring, C., Kohli, P., Jegelka, S., 2018. Batched Large-scale 
Bayesian Optimization in High-dimensional Spaces. In: Storkey, A., 
Perez-Cruz, F. (Eds.), Proceedings of the Twenty-First International 
Conference on Artificial Intelligence and Statistics. Vol. 84. Proceedings 
of Machine Learning Research. 


Wassermann, G., Grewen, J., 2013. Texturen metallischer Werkstoffe. 
Springer-Verlag. 
Weygand, D., Bréchet, Y., Lepinoux, J., Gust, W., 1999. Three-dimensional 


grain growth: a vertex dynamics simulation. Philos Magaz B 79 (5), 
703-716. 


Wicht, D., Schneider, M., Böhlke, T., 2020a. An efficient solution scheme 
for small-strain crystal-elasto-viscoplasticity in a dual framework. 
Computer Methods in Applied Mechanics and Engineering 358, 112611. 


Wicht, D., Schneider, M., Böhlke, T., 2020b. On Quasi-Newton methods 
in FFT-based micromechanics. International Journal for Numerical 
Methods in Engineering online, 1-34. 


Williams, C. K. I., Rasmussen, C. E., 2006. Gaussian Processes for 
Machine Learning. Vol. 2. The MIT press, Cambridge. 


Wilson, J. T., Hutter, F., Deisenroth, M. P., 2018. Maximizing acquisition 
functions for Bayesian optimization. In: 32nd Conference on Neural 
Information Processing Systems (NeurIPS 2018). 


196 


Bibliography 


Wu, B., Vajragupta, N., Lian, J., Hangen, U., Wechsuwanmanee, P., 
Münstermann, S., 2017. Prediction of plasticity and damage initiation 
behaviour of C45E+ N steel by micromechanical modelling. Materials & 
Design 121, 154-166. 


Wulfinghoff, S., Bayerschen, E., Böhlke, T., 2013. A gradient plasticity 
grain boundary yield theory. International Journal of Plasticity 51, 33-46. 


Yang, S., Dirrenberger, J., Monteiro, E., Ranc, N., 2019. Representative 
volume element size determination for viscoplastic properties in poly- 
crystalline materials. International Journal of Solids and Structures 158, 
210-219. 


Zaefferer, S., Wright, S. I., Raabe, D., 2008. Three-dimensional orientation 
microscopy in a focused ion beam-scanning electron microscope: a 
new dimension of microstructure characterization. Metallurgical and 
Materials Transactions A 39 (2), 374-389. 


Zeman, J., Vondřejc, J., Novák, J., Marek, I., 2010. Accelerating a 
FFT-based solver for numerical homogenization of periodic media 
by conjugate gradients. Journal of Computational Physics 229 (21), 
8065-8071. 


Zheng, Q.-S., Zou, Y.-B., 2001a. Orientation Distribution Functions for 
Microstructures of Heterogeneous Materials (I)-Directional Distribution 
Functions and Irreducible Tensors. Applied Mathematics and Mechanics 
22 (8), 865-884. 


Zheng, Q.-S., Zou, Y.-B., 2001b. Orientation Distribution Functions for 
Microstructures of Heterogeneous Materials (I)-Crystal Distribution 
Functions and Irreducible Tensors Restricted by Various Material 
Symmetries. Applied Mathematics and Mechanics 22 (8), 885-903. 


Zhuang, X., Wang, Q., Zhu, H., 2015. A 3D computational homog- 
enization model for porous material and parameters identification. 
Computational Materials Science 96, 536-548. 


197 


Schriftenreihe Kontinuumsmechanik im Maschinenbau 
Karlsruher Institut für Technologie (KIT) 
(ISSN 2192-693X) 


Band 1 


Band 2 


Band 3 


Band 4 


Band 5 


Band 6 


Felix Fritzen 

Microstructural modeling and computational homogeni- 
zation of the physically linear and nonlinear constitutive 
behavior of micro-heterogeneous materials. 

ISBN 978-3-86644-699-1 


Rumena Tsotsova 
Texturbasierte Modellierung anisotroper Fließpotentiale. 
ISBN 978-3-86644-764-6 


Johannes Wippler 

Micromechanical finite element simulations of crack 
propagation in silicon nitride. 

ISBN 978-3-86644-818-6 


Katja Jöchen 

Homogenization of the linear and non-linear mechanical 
behavior of polycrystals. 

ISBN 978-3-86644-971-8 


Stephan Wulfinghoff 

Numerically Efficient Gradient Crystal Plasticity 
with a Grain Boundary Yield Criterion and 
Dislocation-based Work-Hardening. 

ISBN 978-3-7315-0245-6 


Viktor Müller 

Micromechanical modeling of short-fiber 
reinforced composites. 

ISBN 978-3-7315-0454-2 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


Band 7 


Band 8 


Band 9 


Band 10 


Band 11 


Band 12 


Band 13 


Band 14 


Florian Rieger 
Work-hardening of dual-phase steel. 
ISBN 978-3-7315-0513-6 


Vedran Glavas 

Micromechanical Modeling and Simulation 
of Forming Processes. 

ISBN 978-3-7315-0602-7 


Eric Bayerschen 

Single-crystal gradient plasticity with an accumulated 
plastic slip: Theory and applications. 

ISBN 978-3-7315-0606-5 


Bartholomäus Brylka 

Charakterisierung und Modellierung der Steifigkeit von 
langfaserverstarktem Polypropylen. 

ISBN 978-3-7315-0680-5 


Rudolf Neumann 

Two-Scale Thermomechanical Simulation 
of Hot Stamping. 

ISBN 978-3-7315-0714-7 


Mauricio Lobos Fernandez 

Homogenization and materials design of mechanical 
properties of textured materials based on zeroth-, 
first- and second-order bounds of linear behavior. 
ISBN 978-3-7315-0770-3 


Malte Schemmann 

Biaxial Characterization and Mean-field Based Damage 
Modeling of Sheet Molding Compound Composites. 
ISBN 978-3-7315-0818-2 


Jurgen Albiez 

Finite element simulation of dislocation 
based plasticity and diffusion in multiphase 
materials at high temperature. 

ISBN 978-3-7315-0918-9 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


Band 15 Maria Loredana Kehrer 
Thermomechanical Mean-Field Modeling and Experimental 
Characterization of Long Fiber-Reinforced Sheet Molding 
Compound Composites. 
ISBN 978-3-7315-0924-0 


Band 16 Peter Hölz 
A dynamic and statistical analysis of the temperature- and 
fatigue behavior of a race power unit - The effect of differ- 
ent thermodynamic states. 
ISBN 978-3-7315-0988-2 


Band 17 Andreas Prahs 
A Gradient Crystal Plasticity Theory 
Based on an Extended Energy Balance. 
ISBN 978-3-7315-1025-3 


Band 18 Johannes Ruck 
Modeling martensitic phase transformation in 
dual phase steels based on a sharp interface theory. 
ISBN 978-3-7315-1072-7 


Band 19 Hannes Erdle 
Modeling of Dislocation - Grain Boundary Interactions in 
Gradient Crystal Plasticity Theories. 
ISBN 978-3-7315-1196-0 


Band 20 Johannes Görthofer 
Microstructure generation and micromechanical modeling 
of sheet molding compound composites. 
ISBN 978-3-7315-1205-9 


Band 21 Daniel Wicht 
Efficient fast Fourier transform-based solvers for computing 
the thermomechanical behavior of applied materials. 
ISBN 978-3-7315-1220-2 


Band 22 Juliane Lang 
Thermomechanical Modeling and Experimental 
Characterization of Sheet Molding Compound Composites. 
ISBN 978-3-7315-1232-5 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


Band 23 Julian Karl Bauer 
Fiber Orientation Tensors and Mean Field Homogenization: 
Application to Sheet Molding Compound. 
ISBN 978-3-7315-1262-2 


Band 24 Sebastian Gajek 
Deep material networks for efficient scale-bridging in 
thermomechanical simulations of solids. 
ISBN 978-3-7315-1278-3 


Band 25 Jannick Kuhn 
Microstructure modeling and crystal plasticity parameter 
identification for predicting the cyclic mechanical behavior 
of polycrystalline metals. 
ISBN 978-3-7315-1272-1 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


To capture the influence of the microstructure on the cyclic mechanical be- 
havior of polycrystalline metals in a resource-efficient way, computational 
homogenization methods solve equations on computational cells. These reflect 
the essential characteristics of the microstructure. In this work we propose to 
use modern solvers to compute representations of polycrystalline materials in 
terms of Laguerre tessellations with prescribed volume fractions. We assess the 
performance of these solvers in terms of the number of iterations and overall 
run-time needed to compute the Laguerre cells. To assign a crystallographic 
orientation to each cell, we propose an optimization scheme based on tensorial 
texture coefficients, a measure used to quantify the anisotropy of polycrys- 
tals. We investigate the capability of the method to reproduce isotropic and 
anisotropic mechanical behavior with microstructures featuring different grain 
size distributions. Describing the deformation behavior on a microscopic scale 
requires identifying suitable parameters for the material model at hand, i.e., 
single crystal plasticity. As the ground truth we use macroscopic experiments 
and compare them with computational results, leading to an inverse optimiza- 
tion problem. We propose to use a Bayesian optimization strategy, investigate 
its performance as well as its robustness and compare it to other algorithms. 


15-1272-1 


ISSN 2192-693X 
ISBN 978-3-7315-1272-1 
9 "7337314512721 


Gedruckt auf FSC-zertifiziertem Papier 


272 


